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4.0 Executive Summary 

The USM3D Navier-Stokes flow solver contributed heavily to the NASA Constellation Project 
(CxP) as a highly productive computational tool for generating the aerodynamic databases for 
the Ares I and V launch vehicles and Orion launch abort vehicle (LAV). At the peak of that 
effort, USM3D was consuming approximately 25 percent of the total available computer 
resources on the 10,000-processor NASA Advanced Supercomputing Columbia machine at the 
NASA Ames Research Center for these projects. 

Some of the largest database uncertainties were associated with inadequate aerodynamic 
modeling of jet impingement. For example, the aerodynamic modeling uncertainties from the 
abort motor plume impingement on the LAV flight vehicle continually contributed to failure of 
the flight simulations, which elevated the risk to flight safety for launch abort. Similar 
uncertainties were present in the jet impingement from the Ares roll control system (RCS). 

There is a need to improve this deficiency through implementation of more physically correct 
models of the hot-gas jet flows. Efforts to resolve these uncertainties within any of the flow 
solvers used in the CxP have not been successful to date. 

USM3D is currently limited to ideal-gas flows, which are not adequate for modeling the 
chemistry or temperature effects of hot-gas jet flows. This task was initiated to create an 
efficient implementation of multi-species gas and equilibrium chemistry into the USM3D code to 
improve its predictive capabilities for hot jet impingement effects. A finite -rate chemistry model 
was also pursued to improve future predictions of re-entry heating rates of planetary entry 
vehicles. 

The effort was only partially successful as the investigators struggled with converting an existing 
flow solver with tightly embedded ideal-gas assumptions to a real-gas chemistry code. At this 
time, the implementation of the real-gas and multi-species transport equations has been verified 
on simple unit problems, but has not been validated against experimental data. Also, some 
numerical solution stability issues are yet to be resolved. The modules for the equilibrium and 
reacting chemistry were developed, but were not implemented because of the prior difficulties. 

Work is continuing through the support of the Space Launch System (SLS) Aero Task Lead, 

Dr. John Blevins at NASA Marshall Space Flight Center, using the services of subject matter 
experts, including: Dr. Michael Applebaum, Dr. William Eppard, and Mr. Les Hall at CRM 
Solutions, Inc., in Huntsville, Alabama. Assessments are underway to determine the best path 
forward. The near-term goal is to resolve the numerical stability issues with the multi-species 
capability as soon as possible and commence application to Multi-Purpose Crew Vehicle 
(MPCV) launch abort and SLS RCS jet cases. The long-term goal is to implement and test the 
equilibrium and reacting-gas capability. 
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5.0 Assessment Plan 

The goal of this assessment was to implement and validate a simulation capability to handle real- 
gas effects in the USM3D code. This implementation was proposed as a three-step process: 

1 . Modification of governing equations to handle variable ratio of specific heats. 

2 . Addition of species transport equations to model the transport and mixing of multiple 
species in the computational domain. 

3 . Addition of capability for different species to react with each other and produce new 
species. 

A detailed description of the USM3D code modifications, addition of new capabilities, and 
implementation model validation are described in the following sections. 

6.0 Problem Description 

The underlying physics associated with flow around a supersonic/hypersonic vehicle is complex, 
which involves an oblique/bow shock, air dissociation (or ionization), and large heat transfer. 
Accurate predictions of the aerothermodynamics of supersonic/hypersonic vehicles require 
efficient numerical models for resolving chemical reactions and modeling of transport of 
multiple species, and accurate estimation of heat transfer. The modeling of chemical reactions 
and species transport is also important in the numerical simulation of combustion. Challenges 
involved in the simulation of flows with chemical reactions are the stability and efficiency issues 
associated with the numerical methods. The stability of the numerical scheme for the species 
equation is hindered by disparity between the time scales for the mean flow and chemical 
reactions and the stiffness problem of the source terms. Various methods were developed to 
alleviate the stiffness problem and to efficiently solve a large set of species transport equations 
with momentum and energy equations. Since each method has its strength and weakness, the 
USM3D code needed to be analyzed with the assistance of the code developer to determine the 
numerical method most appropriate to efficiently and accurately solve the set of governing 
equations. This report describes the addition of the chemically-reacting flow capability into the 
USM3D code. 
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7.0 Technical Approach 

Section 14.1 provides the nomenclature definitions for equations used in the following sections. 

7.1 Theoretical Formulation 

Governing Equations 

The governing equations for single-species real-gas flow can be expressed as: 

Continuity equation: ^ + V ■ f pV 1 = 0 (1) 


Momentum equation: 
Total energy equation: 
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The reference length scale ( i ref ) was read into the USM3D code through input based on the unit 
of the mesh coordinates. If the air at the sea level (1 atm and 300K) was taken as reference state, 
then the reference temperature ( T ie{ ), pressure ( p re{ ), species mass fractions ( a s ref ), the 

reference gas constant ( i? ref ), reference density ( p Kt ), reference specific heats ratio ( y Kf ), 
reference speed of sound ( a ref ), and reference viscosity ( p ref ) were determined automatically. 

The values of all flow variables and mesh coordinates were non-dimensionalized based on this 
set of reference quantities once they were read into the code. The solution of the flowfield was 
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calculated based on the non-dimensionalized governing equation. The Sutherland law for air 
was used as: 
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The set of non-dimensionalized governing equations was obtained as: 
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For the ease in discretization of these equations for unstructured meshes, the governing equations 
were written in the integral form as: 
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where P is the contravariant velocity component of the fluid particle with respect to the grid and 
is defined as: 


P — im x + vn y + \vn_ — n t 

The thermal equation of state for perfect a gas is: 

t = 7 iSL P ( 14 ) 

pR 

Relations between enthalpy per mole, specific heat per mole at constant pressure, and 
temperature are: 
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Species Continuity Equations 

One of the prerequisites for the chemically-reacting flow modeling is the capability to solve 
species transport equations. This is a result of chemical reactions, where species mass fractions 
can be changed and transported through convection and diffusion processes. The local mass 
fraction of different species with the temperature and pressure govern the local reaction. The 
change of species mass fractions alters the mixture temperature and density at a given pressure 
and enthalpy. Therefore, an accurate modeling of the species transport equation is an important 
step. Appropriate equations of state (i.e., thermal and caloric) need to be employed to account 
for the real-gas effect involving multiple species. 
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The species transport equation can be written as: 


d pa t 
dt 
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(16) 


where p and V are the mean density and velocity vector, D. , « , and ox are diffusion coefficient, 
species mass fraction, and source term of i-th species, respectively. 

Using an eddy viscosity model, the diffusion coefficient can be written as: 
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Integrating over the control volume and applying the divergence theorem, Equation 1 8 reduces 
to: 
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where i = 1,2, NS and NS is the number of species 

For modeling chemical reactions, the species source term can be evaluated based on local species 
mass fractions and thermal state using either equilibrium or finite -rate chemistry. The species 
transport equation can be solved either independently (i.e., segregated approach) or by coupling 
with other governing equations (i.e., coupled approach). For the segregated approach, the 
species equation can be solved with the same time-step size as other governing equations, or with 
a smaller time-step size than other equations due to the disparate time-scale between flow and 
chemical reactions. The segregated approach requires sub-iterations to ensure convergence. The 
coupled approach does not need sub-iterations, but is limited to use the same time-step size for 
all equations, and thus requires large computer memory to store the Jacobian matrix, especially 
when there are a large number of species involved. 
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For the source term of the species transport equation, it can be solved using an explicit or 
implicit approach. Each of these methods has advantages and disadvantages. In the case of 
explicit methods, the source term will be treated explicitly so that it is computationally efficient. 
However, the stability and accuracy are sacrificed due to the stiffness of the source term. In the 
case of implicit schemes, the sources terms will be coupled together and treated implicitly. This 
increases the slow solver stability, but the computational efficiency will be compromised 
because of the need for the solution of the denser matrix system. 


Equilibrium Chemistry Model 

The general form of chemical reactions can be expressed as: 
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bj 
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where v' and v" are the stoichiometric coefficients of species in forward and backward 

J l J l 

reactions K f . and K b . are the forward and backward reaction rates of the f 1 reaction, O, is 
the name of /-th species, and NS is the total number of species. 

The equilibrium chemistry model assumes that the forward and backward reactions are in the 
equilibrium state. In this case, the species mass fractions are determined based on the local 
thermal state (i.e., pressure and mixture energy) and element balance (i.e., chemical elements 
cannot be created or destroyed). Based on these principles, a Gibbs energy minimization method 
can be devised to evaluate the equilibrium species mass fraction ( a " +l ). There are two ways to 

model the equilibrium chemistry with species transport. The first approach is to assume that 
reaction takes place instantly (i.e., infinitesimal time), where the species transport equation is 
solved with zero source term. After the species equation is solved, the equilibrium chemistry is 
used to determine the new species mass fraction based on the local mass fractions and thermal 
state. The second approach is to treat the convection, diffusion, and reaction processes that occur 
simultaneously, and evaluate the source term as: 


/>«-<) 

CO: = 

At 

The chemical equilibrium constant can be used to calculate the reverse rate constant if the 
forward rate constant is known. 


( 21 ) 


Limitations of the equilibrium chemistry include: 

1 . Equilibrium chemistry assumption is valid only at high temperature and high pressure 
when the reaction is extremely fast, 
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2 . Fuel combustion (e.g., heavy hydrocarbon fuel) does not follow the equilibrium 
chemistry, especially at the fuel-rich condition, 

3 . The occurrence of soot and/or nitrogen oxide is controlled by the kinetic chemistry, and 
is not in chemical equilibrium, 

4. Ignition and flame extinctions are not equilibrium processes, and are dependent on 
kinetic rates, and 

5 . The equilibrium chemistry cannot be used to determine the flame front location. 


Finite-rate Chemistry Model 

For a finite -rate chemistry model, which assumes all chemical reactions occur at a finite rate, the 
source term in Equation 19 is calculated based on the reaction processes, their associated kinetic 
rates, local species mass fractions, pressure, temperature, and the presence of a catalyst or 
inhibitor. The source term due to chemical reaction can be expressed as: 
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where M t is the molecular weight of species i, and NR is the total number of reactions. The 

chemical reaction rates in Equation 22 controls the rate of species formation and destruction. For 
example, higher temperature and pressure will result in a faster reaction rate. The reaction rate 
formulated in the empirical Arrhenius form is: 

K f =AT B e x P [-E a /R„T] (23) 


where A, B, E a are empirical constants. As described earlier, there are two ways of calculating 
the species source term. For the explicit method, since it is numerically unstable and inaccurate, 
some researchers employed the penalty function method to minimize the numerical error and 
improve the numerical stability. 

In the method, the time-step size of the species equation is determined by: 
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where (A a) a is the specified allowable maximum species mass fraction change. However, for a 
stiff reaction model, the calculated time-step size can become small, and thus requires a large 
number of sub-iterations. For the implicit method, since all the species involved in the source 
term are coupled and treated implicitly, the Taylor series expansion is employed to linearize the 
source term. Because the solution of a matrix is involved in the implicit scheme, it is less 
computationally efficient, but is more stable and accurate. 

Linearization of the Inviscid Flux 


Representing the conserved variable vector Q as: 
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The contravariant velocity can be represented as: 
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Using Equation 26, the inviscid flux vector can be written as: 
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Equation 27 can be differentiated with respect to Equation 25 to get the inviscid flux Jacobian. 
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Each component of the inviscid flux vector can be represented as A . 


in: 


A = 


i 


*21 


12 


*22 


13 


23 


33 


-42 


Al A 2 A 3 A 


d fi 

SQ k 


which will result 


A , , A,, A,, A,, A, 


*14 


*15 




*24 


*25 


*34 


A 3 i A 32 A 

■a a A,, A 

A41 A 


*35 


A 43 A a „ a 


44 


54 


45 


55 


(28) 


For the linearization of the flux vector, it is assumed that the ratio of specific heats is locally 
constant. Based on this assumption, each term in Equation 28 can be expressed as: 

— n, n r 11 n. 0 


A = 


n x tj)-ud 
n d> - vO 
n_(j)-wd 

PJ 


6 n. 


nX2-y)u + U 
n x v - n , (y - l)w 
n x w-n x (y-\)u 

A. 

V P 


where 


0 — un . + vn„ + wn 


1 7 1 f 2 2 2 

(p = \u + v + w 


) 


n y u - n x {y - l)v 
n,(2-y)v + U 
n„w — n_(y - l)v 


\juO n r 


r A ^ 

v P 


-(y-\)vO 


n : u-n r (y-\)w 
ny - n y (y -\)w 
n , (2 -y)w + U 

A. 

v P 


n x (y- 0 

n y(r~ 0 
«_-(/- 0 


~(y-\)w() yd-n t 


(29) 


(30) 


(31) 


Linearization of Viscous Terms 

The viscous terms involves the gradients of the primitive variables. Therefore, linearization of 
the viscous terms will be simpler if primitive variables are used as the independent variable. To 
linearize the viscous flux with respect to the conserved variables, the Jacobian, based on the 
primitive variables, can be multiplied by the linearization of the primitive variables with respect 
to the conserved variables. 
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The variable vector q selected for this linearization procedure is taken as: 


P 


<h 


a 

u 




a/e, 

V 

= 


= 

e 3 /e, 

w 




e 4 /e, 

p 


<?5 


(■is 


( 32 ) 


Hence, Jacobian resulting from the linearization of the primitive variable with respect to the 
conserved variable can be written as: 


dq 

50 


dq x dq x dq x dq x dq x 


dQ i 

50 2 

503 

50 4 

505 


Sq 2 

dq 2 

dq 2 

5<?2 

dQ x 

50 2 

503 

50 4 

505 

dq 3 

dq 3 

dq 3 

5^3 

5^3 

dQ x 

502 

503 

504 

505 

dq 4 

dq 4 

dq 4 

dq 4 

dq 4 

50, 

502 

503 

504 

505 

dq 5 

dq 5 

dq 5 

5^5 

5^5 

50, 

502 

503 

504 

505 


( 33 ) 


Using the relation between primitive and conserve variable, the above Jacobian can be reduced 
to: 


dq 

dQ 


1 0 

-u/ p Up 

-v/ p 0 

-wl p 0 

§* -O'- D« 

dp 


0 0 0 

0 0 0 

1 Ip 0 0 

0 Up 0 

-(y-l)v -{y -\)w y-l 


( 34 ) 


The details of this derivation are provided in Appendix A. The linearization of pressure for real- 
gas is provided in Appendix B. 


NESC Request No.: 07-037-1 



NASA Engineering and Safety Center 
Technical Assessment Report 

Document #: 

NESC-RP- 

07-037 

Version: 

1.0 

Title: 

Reacting Multi-Species Gas Capability for USM3D Flow 

Solver 

Page #: 

19 of 78 


From Equation 14, Equation 35 can be obtained as: 


dT = -—dp + ^^d p 
P P R 


Hence, 

dT 

dq 


— 0 0 0 

P P R 


(35) 


(36) 


Combining Equations 34 and 36 produces: 


dT 

dQ 


8T 8T 8T dT dT 


dQi oQ 2 dQi oQa d Q 5 

1 


— 0 0 0 ^ 
P P R 


dT dq 




dq dQ 




0 

0 

0 

0 

Up 

0 

0 

0 

0 

Up 

0 

0 

0 

0 

Up 

0 

(y-l)u 


-(y-l)w 

r- 1 


-u/ p 
-vl p 

-wl p 
dp 

dp 

-T + 7 refdp -r t Ar-Vu -r re[ (r-i)v -r t Ar-Vw twO'-i) 


p pR dp 


P R 


P R 


P R 


P R 


-t + r r ef (r~i)v 2 y ref (y i) w - y ,A y - 1) v -YtAy- l ) w Y xAY ~ 1) 


p pR 2 


pR 


pR 


pR 


P R 


(37) 


The spatial derivatives any primitive variable (ft at the interface between two adjacent cells can be 
calculated as: 


(j) x = ~ = E v ; T x = ^ 

dx n x Ax + n y Ay + n z Az 

<f> =^-=r a</>; r y = ^ 

dy ' ' n x Ax + n y Ay + n_Az 

(j) = = r A(j ) ; ^ 

dz ‘ " n r Ax + n Ay + n.Az 

x y ^ z 


( 38 ) 
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and A(/) = (f) R — (/> L where (f> R and <f) L are the values of <j) at the center of right and left cells, as 
shown in Figure 7.1-1. 



Figure 7.1-1. Nomenclature of a Face and it's Neighbors 

Hence, the spatial derivatives of a dependent variable at the cell face can be linearized by the 
conservative Q as: 


M. =r a(A0 _ r 

3Q x 8Q \ dQ , 


(39) 
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Therefore 


SQ r 

d^_ 

dQ L 


=r 


K d Q R J 

r 


(40) 


= -r 


8<t>L 
k q Qlj 


(41) 


Equation 13 can be differentiated with respect to Equation 25 to get the viscous flux Jacobian. 


Each component of the viscous flux vector can be represented as B jk = 


dJl 

SQ k 


and this will result 


in: 


B = 


1 ^12 ^13 ^14 B 15 

^2i B 22 B 23 B 24 B 25 

B 3i B 32 B 33 B 34 B 35 

B41 B A2 B 43 B 44 B 45 

B51 B 52 B 53 B 54 B 55 


(42) 


By substituting Equations 34, 37, and 39, each components of the viscous flux Jacobian can be 
expressed as: 


B\ 1 — B X2 — B a — B l4 —B l5 — 0 


o _1IL 


s/; 


3r„ 


8r„ 


-+n 


+n„ 


dr. 


Bjj = Xqc), 


B « = 


3/; 


3(a), ' 3(a), : 3(a), 


3r„ 


3r. 


dr. 


-+n 


+n_ 


3(a), y 3(a), z 3(a),- 




3r„ 


3r. 


3(a), I x 3(a), y d{Q c )j z d{Q c ) 

8fs _ d ( u .ffi) , g(v/;) . 3(1 vf 4 ) 


(43) 

(44) 

(45) 

(46) 


Bsj 3(a),. 3(a),. + 3 (a)! + 3(aV +/r 


8T 


8T 


-+11 


+n_ 


8T 


s(Q c ), ’s(Qc), ' d (Qc), J 


(47) 


= “/ fl 2y +V / fl 3,'+ W / fl 4y+* : 


5T 


-+W 


sr 


+n 


57 


0(2c), ’0(2c), '0(2c)J 
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where C represents either cell on the left or the cell on the right-side of the face and j represents 
the column number. 

7.2 Numerical Implementation 

An overall flow diagram of the USM3D code main time loop with the current modifications to 
add real-gas capability is given below. The sections highlighted in bold letters represent new 
code additions. The USM3D code was modified to add variable ratio of specific heats. The 
modifications added to the different code modules are described in this section. 

DO 25 ntsUN = 1, ntstep 
Move grid 
Qnode 
UBC2 

DO 10 nts = f niter 
UCTime 
Limiter 

Turbulence model -> Update viscosity 

Calculate species source terms based on equilibrium or finite chemistry 
Solve species transport equations 
Save (pa )" +1 

Viscous flux 

Update (Implicit scheme ) 

Calculate RHS 
Calculate LHS 
Calculate AQ 
Update state variables 

Update species mass fraction a"' = x 

P 

Calculate new species mass fraction based on equilibrium chemistry 
model if the effect is not included in the species source term 
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Update temperature and pressure using new species mass fractions and 
state variables 

Qnode 

UBC2 

Force 


10 Continue 
25 Continue 

Roe’s Approximate Riemann Solver 

Based on Roe’s approximate Riemann solver, the numerical flux that crosses a cell-face can be 
written as: 


F LR 2 


((ffoj+ffe,)- a|(& - ej) 


(48) 


where L and R represent the cells on the left and right sides of the face respectively, Q is the 
conserved variable vector, and A is the Roe averaged matrix. The flux vector F{Q) is defined 


as: 


F(Q) = 


pU 

puU + pn x 
pvU + pn y 
pwU + pn z 
U(pE + p)_ 


(49) 


where E is the total energy per unit mass and U is the contravariant velocity. The variables used 
in the Roe averaged matrix are evaluated as: 


P — • J PlPr 

- \Ul +u r^Pr/Pl) 

i V Pr ( Pl 

- Vl+VrtIPrI Pl) 

v = , — - 

1 + sIPr/Pl 


(50) 

(51) 

(52) 
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- ( w l + w r^Pr / Pl) 
w = — - 

! + a IPr/Pl 

T _ (/ 7 0L hpR V Pr ! Pl) 
\ + ^Pr! Pl 


(53) 


(54) 


The diffusion term (i.e., last term in Equation 48) can be split it into three components as: 

\a\(Q r - Gl ) = | at; | + 1 af 2 | + 1 af, | 


(55) 


where, AF, = \u\ 


a J 


1 

u 

V 

w 

Hu 2 +v 2 + w 2 ) 
2 V ; 


P 


0 

Au-n x AU 
Av - n y AU 
Aw-n_AU 

uAu + v Av + w Aw — UAU 


1 



\U ±a\ 


r Ap±pa AU^ 

v 2d" 2 , 


u + n x a 
v + n y a 
w ±n,a 


h Q ±U a 


(56) 


(57) 


where A represents the change in the variable value between the cell on the right and cell on the 
left. 

The different components (i.e., mass, momentum, and energy) of the flux vector can be written: 


(pU) L + (pU\ P\( A A 

J \ /-\ /-\ / — 2 

2 2 V a J 


— \U + a\ 


f \ 

A p+ paAU 


4 a 


— 2 


-U-a\ 


c ^ 

Ap-paAU 




4a 


—2 


(58) 
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fl = 


C puu\ + (puU) R + (p, +p R )n x | tv 7 1 


— \U + a\ 


f A 

A p+ paAU 


4 a 2 


2 2 

(u + n x a)- 1 U - a 


u\Ap-^Y + p(Au-n x AU) 


y ■\ 

Ap-paAU 


4a 2 


( u - n x a ) 


A = 


( P vU )l + (p vU )r + (Pl+ Pr K M 
2 2 

y A 

\YJ _| Ap+paAU f 1 = 

-|t/ + a| — ^ (v + n v a\-\U - 


y y 


v 

V v 


A/? - + p(Av - n y AU)\ 


4a 2 


f A 

Ap-paAU 

4 ¥ 


(v - « v a) 


f* = 


(pwU\ - (pwU) R + (p L + p R )n z |t/| 
2 2 

(w + n : a^-\U - 


w 


-P- 


■ a\ 


^ A p + paAU^ 


4 a 


a\ 


Ap— + p(Aw — n z AU) 

v a ' J 

^ Ap- pa AU^ 


4 a 


— 2 


(w-n z a) 


h=~ 


(pU(E + p)\ + (pU(E + p)) 

^PP}_Er 2 , -2 , — 2M — 


£/ 


Ap - 1 -^(u 2 + v 2 + w 2 )l-p(w Aw + vAv + wAw-U At/) 


£/ + G 


^ Ap + paAU^ 


r7-\ It T A P~P aAU 


/ 5 = 


2a 2 

(pUh 0 ) L +(pUh 0 ) R \U 


(h 0 + Ua)-\U 


A P~^r 
a y 


— icy — a 
Ap 


v 2g 2 

A F 2 |C7| _ 


|^ 0 + £/a) 


— t/ + G 


y X 

Ap + paAU 

4 a 2 


2 2 

y 

+ £/g)— | t/- g| 


p(gAg + vAv + wAw — UAU^ 
Ap-paAU ^ 


4g 2 


(A 0 -t/a) 


(59) 


(60) 


(61) 


(62) 


(63) 


For flows with multiple species, the ratio of specific heats used in the Roe averaged matrix is 
evaluated in different ways by different authors. Two of the commonly used methods are 
described in the following sections. 
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Methods used by Luke [ref. 2] 

The ratio of specific heats is evaluated as: 

- i R 

/ = 1 + (64) 


where, 


_ NS 
s=l 

(65) 

NS 

C*=YccC* , and 

V / J X vs ’ 

.9=1 

(66) 

c;=-Vi;x^ 

(67) 

(Ps,L ^ V Pr! Pl ) 

i+Va/a 

(68) 

t _(Tl + T r ^Pr/Pl) 

\+4p r !pl 

(69) 

(f*L +e *W Pr/Pl) 
1+ ^Pr/Pl 

(70) 

The speed of sound is calculated as: 


« 2 = (r - 4 /7 0 - - (^ 2 + )+ - J^a s e s 

V Z S=1 J 

(71) 

Method used by Molvik and Merkle [ref. 3] 


Molvik and Merkle used the following approach to calculate the speed of sound and the ratio of 
specific heats. The initial estimation for the speed of sound, ratio of specific heats, and the 
chemical source term are defined as: 

a, +a R 
a = L K 
2 

(72) 

7 =r L + r, 

1 o 

(73) 
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D D sL + D sR 


These values are refined to meet Roe’s criteria as: 


(74) 


- = D s -a Ap s Sp! x 


r = 


l-A pSp! x 

Y - 1 


l-A p5p! x 


+ 1 


(75) 

(76) 


where A represents the change in the value between right and left states, and the other variables 
appearing in Equation 76 are calculated as: 

NS 

5p = A p + Y,D S A p s -{/- l\A(ph)-Ap) (77) 

5=1 


NS . 

X = A A )‘ + (Sp) 2 (78) 

5=1 

The Roe averaged speed (a) of sound is calculated as: 

_ NS 

a 2 =(y- 1 )h - YjP s D s (79) 

5=1 

Methodology Implemented in USM3D Code 

Currently, in the USM3D code, the Roe averaged speed of sound is calculated based on constant 
specific heat assumption expressed as: 

a 2 = /RT =y(y-l)C v T = -^(u 2 +v 2 +w 2 )j (80) 

This has been modified to calculate the speed of sound based on variable ratio of specific heats. 
The Roe averaged density, velocity components, and total energy were calculated using 
Equations 50 through 54. The other Roe averaged correlations used for different flow variables 
are listed as: 


a „ 


(^ s ,L+«s,R-JPr/Pl) 

1 +^Pr/Pl 

E ——(u‘ 


e =E — -\u 2 +v° + w’ ) = 

5=1 


( 81 ) 

( 82 ) 
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Since de = C v dT , the Newton-Raphson method was used to solve for T(=T'" ) as: 

— n + 1 —n 

n + 1 jTfl ^ ^ 

C" 

v 

— n + 1 —n 

e e 
Y n+ i y n ^ 

C " /r 


S = S^ = S_i = S_ 1 = y„- ( C A 

R D D D s 


R 


R 


R. 


R. 


-1 


e h-RT h - h ° - ^ /?° - 

- — = — = = -T= T=Yn T 

R R R R„ tT R 


a.lM. 


n „ = 


J AS 


Z(«./aO 

5=1 

(c p) j R u — {ciyT + ci^T + ci 3 + ci + ci^T -\-ci^T -\-q,-jT 

/z v j R u =(-a,r +<7 2 1 nr +a^T + y ct^l' -\--^cl 5 T~ + u ( T + 4 u^T +£t^ 
Once the Roe-averaged temperature is obtained, a , can be calculated based on: 
= A SfRT 


a 

where 


AS 




«, = — 

M. 


C./i? 

' (C //?)-! 


(83) 

(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 

(90) 

(91) 

(92) 

(93) 
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Modifications to Boundary Conditions 

Boundary Condition 0: Freestream (Supersonic Inflow) 

• Specify and fix density, velocity components, species mass fraction (i.e., mixture gas 
constant), and pressure. 

• Calculate mixture temperature and Gamma (i.e., ratio of specific heats) based on 
pressure, mixture density, and mixture gas constant using the thermal equation of state. 

• Calculate total energy from species mass fractions and mixture temperature, and velocity 
based on the caloric equation of state. 

Boundary Condition 1: Tangent Flow (Reflection Plane) 

• Set contravariant velocity component to zero. 

• Set pressure, mixture density, and species mass fractions at the boundary to be the same 
as the values of the cell values. 

• Calculate mixture temperature and Gamma based on pressure, mixture density, and 
mixture gas constant using the thermal equation of state. 

• Calculate total energy from species mass fractions and mixture temperature, and velocity 
based on the caloric equation of state. 

Boundary Condition 2: Full Extrapolation (Supersonic Outflow) 

• Extrapolate velocity components, pressure and mixture density at the boundary from the 
cell values. 

• Set species mass fractions at the boundary to be the same as the values of the cell inside. 

• Calculate mixture temperature and Gamma at the boundary based on pressure, mixture 
density, and mixture gas constant using the thermal equation of state. 

• Calculate total energy at the boundary from species mass fractions and mixture 
temperature, and velocity based on the caloric equation of state. 

Boundary Condition 3: Characteristic Inflow/Outflow (Subsonic Outer Boundaries) 

Appendix C provides a detailed description. 

• Assumptions 

o Species mass fractions at the boundary are the same as species mass fractions inside if 
the contravariant velocity is positive (i.e., outflow), 
o Species mass fractions at the boundary are the same as freestream species mass 
fractions if the contravariant velocity is negative (i.e., inflow). 

• Use Riemann invariants to calculate the boundary contravariant velocity and the speed of 
sound. Use the speed of sound to calculate the temperature. 

The Riemann invariants corresponding to the outgoing and incoming characteristics are given by: 

9T = U + — , and 9T = U - — 
y-l y — 1 
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where S J? is corresponding to the eigenvalue, U+a, and S H is corresponding to the eigenvalue 
U-a. These Riemann invariants at the boundary can be written as: 


9?+ 

* 1 inside 

if 

U + a> 0 

(94) 

911 „ 

. nee s tie am 

if 

o 

V 

a 

+ 

9T .. 

ins ide 

if 

U-a> 0 

(95) 

^fieestieam 

if 

U-a< 0 


The boundary contravariant velocity ({//,) and speed of sound ( a h ) can be calculated by solving 
Equations 94 and 95. The boundary species mass fractions ( (cc s ) b ) will be the same as those of 
the cell at the upwind side. Since: 

«/, = \ YbRhTb ; R h = f(.oc sb ) ; and y h = f(T h , ( a s ) b ) 

V Y*f 

An iterative procedure was used to determine the boundary temperature and the ratio of specific 
heats based on the boundary speed of sound and species mass fractions. With the isentropic 
assumption between the upwind cell and the boundary, the boundary mixture density and 
pressure can be calculated based on the temperature (T c ) and mixture density ( p c ) of the upwind 
cell as: 


i 



Total energy at the boundary was calculated from boundary temperature, species mass fraction, 
and velocity components based on the caloric equation of state. 

Boundary Condition 4: No-Slip (Viscous Surfaces) 

Appendix D provides a detailed description. 

• Stationary wall: set velocity components at the wall boundary to be zero and velocity 
components at the ghost point to have the opposite direction to those at the neighboring 
interior cell. 

• Moving wall: set velocity components at the wall boundary to be the prescribed speed of 
the moving wall and extrapolate velocity components at the ghost point from the 
neighboring interior cell and wall boundary. 

• Extrapolate pressure at the wall boundary from neighboring interior cells, and extrapolate 
pressure at the ghost point from the neighboring interior cell and wall boundary. 

• Set the species mass fractions to be the same as those of the neighboring interior cells. 
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• Adiabatic: in the current implementation, temperature at the wall boundary is set to the 
same as its neighboring interior cell. 

• Isothermal: temperature at the wall boundary remains unchanged. 

• Calculate mixture density: total energy at the wall boundary is calculated from 
temperature, pressure, species mass fraction, and velocity components based on the 
thermal and caloric equation of state. 

Boundary Condition 5: Tangent Flow (Inviscid Aerodynamic Surface) 

• Same as boundary condition 1. 

Boundary Condition 44: Special no-slip (Blunt base) 

• Same as boundary condition 4. 

• If wall function, then it is different than boundary condition 4. 

8.0 Verification of Multi-Species Capability 

The implementation of the numerical models from Section 7.0 was verified using hot-gas flow 
conditions though a round convergent-divergent nozzle exiting into ambient air. The solutions 
presented in the following sections will only serve as qualitative verification of the 
implementation of the real-gas energy and multi-species transport equations, and are not 
validated against experimental data. 

The configuration consists of a plenum chamber, which is represented as the small straight 
section on the left-side of the geometry shown in Figure 8.0-1, followed by the convergent 
divergent nozzle. A computational grid of 2,706,755 tetrahedral cells was constructed for one 
quadrant of the round nozzle configuration. Figure 8.0-1 depicts the grid clustering on one of the 
quadrant boundaries. The hot jet flow exits the nozzle at nominally Mach 2 into near-quiescent 
ambient air. The freestream Mach number is set to 0.05 for numerical purposes. The plenum is 
set to total pressure, p tot , of 1 15 psi and total temperature, T tot , of 2009 degR to produce a hot jet 
plume. The chemical composition of the plenum gas is varied to simulate different flow 
scenarios. The inflow boundary condition to the plenum is specified using Engine boundary 
condition in USM3D code. A description of the USM3D input parameters is presented in 
Appendix E. 
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M e*it = 2 

Figure 8.0-1. Triangulation on Boundary Plane of Tetrahedral Quadrant Grid for Nozzle Test Case 


8.1 Comparison of Ideal and Real Gas Simulations 

This section presents the comparison of the simulation results from the ideal- and real-gas 
versions of the USM3D code. The results confirm the expected variation of the ratio of specific 
heats with temperature from the modified energy equation. The multi-species transport 
equations are not solved in this example. In the following comparison, the real-gas is taken as 
air with a mixture of 73.5 percent nitrogen (N 2 ) and 26.5 percent oxygen (O 2 ). This mass 
fraction is assumed to be the same for the entire flowfield. In these simulations, the Courant- 
Friedrichs-Lewy (CFL) number is ramped from 1 to 10 in 800 iterations. The ideal-gas 
simulation was run for 9,000 iterations. By using 9,000 iterations for real-gas simulation, the 
flow was not fully developed and contour deviations for temperature and Mach number near the 
axis of rotation between the ideal- and real-gas simulations were noted. Therefore, the real-gas 
version of the USM3D code was run for 30,000 iterations. The comparison of the pressure, 
temperature, Mach number, and ratio of specific heats is given in Figure 8.1-1. This figure 
shows that contours of pressure, temperature, and Mach number are similar for both ideal- and 
real-air simulations. The variation of ratio of specific heats for the real-gas simulation can be 
seen in Figure 8.1 -1(d). As expected, the ratio of specific heats is lower in the settling chamber 
where the temperature is high, and defuses toward the ambient gamma ratio of 1 .4 as the 
temperature of the jet plume decreases toward the ambient condition. 
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(a) Pressure 



Temperature 



(b) Temperature 
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Mach 



Mach 



(c) Mach Number 


Gamma 



Gamma 



(d) Ratio of Specific Heats 

Figure 8.1-1. Comparison of Results from Ideal and Real Air Simulations 
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8.2 Comparison of Flux Schemes with N 2 Exhaust Gas 

This section presents the comparison of the simulation results using Roe and van Leer flux 
schemes from the real-gas version of USM3D code. In the following comparison, the ambient 
gas is taken as a mixture of 73.5 percent N? and 26.5 percent 0 2 . The plenum/nozzle gas is 
prescribed as 100 percent N 2 . The species transport was solved to propagate the species into the 
domain, along with the energy equation that modifies the ratio of specific heat with temperature. 
Solutions were started with first-order spatial differencing and switched to second-order. A 
comparison of density, pressure, temperature, velocity magnitudes, Mach number, and species 
mass fractions are shown in Figure 8.2-1. The species mass fraction plots show the N 2 mass 
fraction inside the settling chamber is 1 .0, and 0 2 mass fraction is zero. It is observed that the 
distributions of flow variables computed by the Roe scheme and the van Leer scheme are in 
close agreement, thereby increasing confidence in the implementation of the species transport 
equations. An expected non-reacting transport of the N 2 species from the engine exhaust into the 
ambient air was confirmed. 


Air (N2 ,Q2) 
N2 > — 


Roe Flux Difference Splitting Scheme 



Density 

2.50 
1.95 
1.40 
0.85 
0 30 


van Leer flux scheme 



(a) Density 
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Pressure 




(b) Pressure 




(c) Temperature 
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(d) Velocity Magnitude 


Mach Number 


2.10 
1 .57 
1.05 
0.52 
0.00 




(e) Mach Number 
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Gamma 




(f) Ratio of Specific Heats 


Mass fraction 




(g) Species Mass fraction: N 2 
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Mass fraction 




(h) Species Mass fraction: O 2 

Figure 8.2-1. Results from the Simulations using N 2 as the Exhaust Gas 


8.3 Results from Simulations using Hydrogen (H 2 ) as the Exhaust Gas 

This section presents the results from the simulation of the nozzle jet with IT as the exhaust gas. 
With Hi exhaust gas, the equations are extremely stiff and required more effort to run the 
simulation. First, the solution was run with first-order differencing with a maximum CFL 
number of 5. The CFL ramping was set from 1 to 5 over the first 800 iterations. Second, the 
flow variables inside the plenum were initialized to the hot-gas Ho conditions to facilitate 
solution start up. While USM3D has a general capability for modifying prescribed cells, the gas 
chemistry extension was not added to that feature. For the current testing, the x-coordinate of the 
cell centroids was used to flag the cells in the plenum and initialize the variables. The number of 
iterations was set as 30,000, while maintaining the maximum value of CFL number as 5. The 
distribution of density, pressure, temperature, ratio of specific heats, and species mass fractions 
after 30,000 iterations are shown in Figure 8.3-1. This figure shows that the flow is not fully 
developed at the axis of symmetry due to the cell size. Additional iterations are needed to 
converge the solution and clean up the slight turn in the contour lines near the axis. 
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density 
.1 .440 
1.087 
10. 735 
JO. 383 
Iq.031 



Pressure 

1 5.745 
4.460 
3.175 
1.890 
0.605 



(b) Pressure 


Temperature 

D 3.769 
3.059 
2.349 
1.639 
0.930 



(c) Temperature 
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Mach Number 

1 2.088 

1.566 
1 .044 
0.522 
0.000 



(d) Mach Number 


Gamma 

D 1.403 

1.394 

1.385 

1.376 

1.367 



(e) Ratio of Specific Heats 


H2 

1.000 

0.750 

0.500 

0.250 

0.000 



(f) H 2 Mass Fraction 
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N2 

8 0.767 
0.575 
0.384 
0.192 
0.000 



(g) N 2 Mass Fraction 

02 

! 0.233 
0.175 
0.116 
0.058 
0.000 



(h) 0 2 Mass Fraction 

Figure 8.3-1. Results from the Simulations using H 2 as the Exhaust Gas 


9.0 Findings, Observations, and NESC Recommendations 

9.1 Findings 

The following findings were identified: 

F-l. The implemented real-gas energy equation solver enables the USM3D code to account 
for the effect of variable specific heats. 

Solutions require a smaller numerical stability parameter (CFL) to maintain numerical 
stability, which prolongs run times than its ideal-gas counterpart. 
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F-2. The van Leer and Roe flux schemes were modified in the USM3D code for real-gas 
energy equation solver, with the van Leer flux scheme being the most stable. 

F-3. The species transport equation enables USM3D code to model the multiple species effect. 

A substantially reduced CFL number is needed to maintain numerical stability. The 
extent of CFL number reduction depends on the species system. 

F-4. The numerical modules for solving chemical reactions were developed, but not 

implemented and tested because of the unresolved numerical stability issues for some 
multi-species mixing problems. 

F-5. Continuation of this work is being supported by the SLS Program Aero Task Lead 

through the services of Dr. Michael Applebaum, Dr. William Eppard, and Mr. Les Hall 
of CRM Solutions, Inc. 

9.2 Observations 

The following observations were identified by the Investigators: 

0-1. The original USM3D flow solver was tightly embedded with the ideal-gas assumption 
making it difficult to locate every occurrence during the conversion to a real-gas model, 
thereby increasing the chance for inconsistencies. 

0-2. Redundancy in the numerical flux scheme coding led to difficulties in implementing a 
consistent modification for real-gas equations. 

0-3. The USM3D code has several specialized boundary conditions that presented a challenge 
to adapt to real-gas formulation. 

0-4. USM3D uses a one-dimensional array data structure throughout that requires one to have 
an intimate knowledge of the entire code for making modifications. A more modular 
approach with locally defined arrays would be desirable. 

9.3 NESC Recommendations 

The following NESC recommendations are directed to the NASA Technical Fellow for 

Aerosciences and the Aero thermodynamics Community: 

R-l. Seek recommendation from CRM Solutions, Inc., on best approach for maturing a robust 
multi-species and gas chemistry capability in USM3D. 

(F-l through F-5, 0-1 through 0-4) 

R-2. Monitor the SLS Program support of the USM3D code modification for application to the 
MPCV Program crew module launch abort analysis. (F-5) 
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10.0 Alternate Viewpoints 

There were no alternate viewpoints identified during the course of this assessment by the NESC 
team or the NRB quorum. 

11.0 Other Deliverables 

There were no other deliverables during the course of this assessment. 

12.0 Lessons Learned 


No applicable lessons learned were identified for entry into the NASA Lessons Learned 
Information System. 

13.0 Definition of Terms 


Corrective Actions Changes to design processes, work instructions, workmanship practices, 
training, inspections, tests, procedures, specifications, drawings, tools, 
equipment, facilities, resources, or material that result in preventing, 
minimizing, or limiting the potential for recurrence of a problem. 

binding A conclusion based on facts established by the investigating authority. 


Lessons Learned 


Observation 


Problem 


Knowledge or understanding gained by experience. The experience may 
be positive, as in a successful test or mission, or negative, as in a mishap 
or failure. A lesson must be significant in that it has real or assumed 
impact on operations; valid in that it is factually and technically correct; 
and applicable in that it identifies a specific design, process, or decision 
that reduces or limits the potential for failures and mishaps, or reinforces a 
positive result. 

A factor, event, or circumstance identified during the assessment that did 
not contribute to the problem, but if left uncorrected has the potential to 
cause a mishap, injury, or increase the severity should a mishap occur. 
Alternatively, an observation could be a positive acknowledgement of a 
Center/Program/Project/Organization’s operational structure, tools, and/or 
support provided. 

The subject of the independent technical assessment. 


Proximate Cause The event(s) that occurred, including any condition(s) that existed 
immediately before the undesired outcome, directly resulted in its 
occurrence and, if eliminated or modified, would have prevented the 
undesired outcome. 


Recommendation An action identified by the NESC to correct a root cause or deficiency 

identified during the investigation. The recommendations may be used by 
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the responsible Center/Program/Project/Organization in the preparation of 
a corrective action plan. 

Root Cause One of multiple factors (events, conditions, or organizational factors) that 

contributed to or created the proximate cause and subsequent undesired 
outcome and, if eliminated or modified, would have prevented the 
undesired outcome. Typically, multiple root causes contribute to an 
undesired outcome. 

14.0 Acronyms List 

CFL Courant-Friedrichs-Lewy 

CxP Constellation Project 

FT Hydrogen 

LaRC Langley Research Center 

LAV Launch Abort Vehicle 

MPCV Multi-Purpose Crew Vehicle 

Ni Nitrogen 

NESC NASA Engineering and Safety Center 

NRB NESC Review Board 

0 2 Oxygen 

RCS Roll Control System 

SLS Space Launch System 

14.1 Nomenclature 

A Inviscid flux Jacobian or empirical constant in the Arrhenius form of chemical 

reaction rate 

Ajk Components of the inviscid flux Jacobian 

B Viscous flux Jacobian or empirical constant in the Arrhenius form of chemical 

reaction rate 

Bjk Components of the viscous flux Jacobian 

C p Specific heat at constant pressure 

C v Specific heat at constant volume 

dV Elemental volume 

dS Elemental surface area 

E Total energy per unit mass 

e Internal energy per unit mass 

E a Empirical constant in the Arrhenius form of chemical reaction rate 

F' Inviscid flux vector 

Components of inviscid flux vector 

F v Viscous flux vector 
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rv rv 

J 1 ’ J 2 ’ ’ * 

-,/ 5 v Components of viscous flux vector 


h 

Enthalpy per unit mass 


h 0 

Total enthalpy per unit mass 


K »j 

Backward reaction rate of the j th reaction 


K fj 

Forward reaction rate of the j th reaction 


M w 

Molecular mass 


M 

Mach number 


h 

Outward pointing unit normal 


NR 

Number of reactions 


NS 

Number of species 


n, 

Contravariant velocity component of the grid speed and is defined as 

n , = x,n x + y t n y + z,n z 


Components of the unit vector h 


p 

Pressure 


Pr 

Prandtl number 


Q 

Conserved variable vector 


Qu Qi, ■■ 

, £?5 Components of the conserved variable vector 


R 

Mixture gas constant 


Re 

Reynolds number 


R„ 

Universal gas constant 


S 

Sutherland constant 


T 

Temperature 


t 

Time 


U, V, w 

Components of velocity in v, y, and z directions, respectively 


Xt, y t , Zt 

Components of grids speed in x, y, and z directions, respectively 


Greek Letters 


a 

Species mass fraction defined as the ratio of species density to total density 

P 

Contravariant velocity component of the fluid particle with respect to the grid and is 
defined as (5 = 0 — n t 

r„r,.r s 

Geometric factors 


K 

C n 

Coefficient of thermal conductivity and is defined as k = — 


7 

Ratio of specific heats 


P 

Density 


P 

Dynamic viscosity 


CO 

Source term due to chemical reaction 
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6 Contravariant velocity component of the fluid particle and is defined as 

6 — un x + vn y + wn ? 

v Stoichiometric coefficient 

t xx , r xy , t xz , r yy , r yz , t zz Component of the rate of shear stress tensor 
Subscripts 

b Backward reaction 

/ Forward reaction 

ref Reference conditions 

s Species index 

x, y, z Components of a vector along x, y, and z directions, respectively 
oo Freestream conditions 

Over 

A Dimensional quantities 

— > Vectors 

=r Tensors 
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Appendix A. Derivation of Inviscid and Viscous Jacobians 


Representing the conserved variable vector Q as: 


P 


<2i 

pu 


a 

pv 

= 

a 

pw> 


a 

J P E _ 


a. 


Q = 


The contravariant velocity can be represented as: 

0 = ^-(Q 2 n x + Q 3 n y +Q 4 n z ) 

Using Equation (5), the inviscid flux vector can be written as: 

~n t Qi+n x Q 2 + n y Q 3 +n y Q 4 
%(- n,Q\ + n x Q 2 + n Q . 3 + n y Q 4 )f n x p 


F‘ -n = 


// 
f 2 
fi 

n 

fi 


Q\ 

n,Q x + n x Q 2 + n y Q 3 + n y Q 4 )¥n y p 

7f(- n ,Qi + n xQi + * y Q 3 + «v04 h»,p 
£ l\ 

^-{~n t Q x +n x Q 2 + n y Q 3 +n y Q 4 y^-(n x Q 2 +n y Q 3 +n y Q 4 ) 
\l\ \2\ 


(A.1) 


(A.2) 


(A. 3) 


Equation (A. 3) can be differentiated with respect to A.l to get the inviscid flux Jacobian. Each 


component of the inviscid flux vector can be represented as A , 


dfj 

dQ k 


and this will result in: 


A = 


An A 12 A 13 A 14 A 15 

A A A A A 

^21 ^22 ^23 ^24 ^25 

A 31 A 32 A 33 A 34 A 35 

4 4 4 4 4 

^41 ^42 ^43 ^44 ^45 

.A 51 A 52 A 53 A 54 A 55 


(A.4) 
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A = 

-n, n x n y 

n x (j)-u6 n x (2 — y)u + U n v u — n x (j — \)v 

n v </) — v9 n x v — n(y — l)u n v (2 — yjv + U 

n,(j)-w9 n x w — n_(y — \)u n w — n(y — lV 

/ \ ( \ ( N 

Q n x —~(j> -(y-l)u9 n — - cj> -(y- 1>6> 77. 
.V P ) V P ) v J V P ) v 7 

77, 0 

n z ii-n x {y-\)w n x (y- l) 

ny-n y (y-\)w n y (y- 1) 

77 _ (2 - ;/)vv + [/ 77 _ (7 - 1) 

— - (j> —(y — \ )wO yO-n 
VP J 


and 


, dfi A dfl 
A u = — =-n t ; A 2 = ^~ 
dQ dQ 2 


A d fl A 3 fl 

n; A n = — - = n v ; A 4 = — L 
dQ , " 14 dQ t 


n ; A 15 = ^- = 0 

dQ, 


An 


An 


A 3 


dfi d 


dQi dQi 


q; q 7 q 

n t Q 2 + n x ^ + n 

Qx Q x 


+ n. 


Qi Qa 
Q x 


+ n x p 


Q 


= - — 0 + ft,.— 


dp 




Qx 


dp 


- P+n x {2 - y )% 


30 2 


Qx 


dfl 0, , ,X03 

— = ft v — -ft^-l) — 

dQ, y Qx Qx 


dfl 


dQ, 


02 

Qx 


A 24 = ^ = n zZ f-n x (y-V)- 4 


04 

0i 


25 -^ = n x (y-l) 


— 


As x 


A 2 


dQ, 

dfl = _d_ 

dQx dQi 

dfi 


- n ,Q 3 + n 


0203 


03 


+ n v —^ + n, 

Qx y Qx ‘ 


03 04 
01 


+ n x p 


= -9i g+n l p 

Qx dp 


dQ 2 


03 , n 02 


0, 


0, 


33 - ff = /?+«, (2 -rA 

01 


— 


A 4 


As 


303 

M 

304 
df[_ 

305 


03 , n 04 


0i 

= « y (r-i) 
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Al 


df : d 


dQ, 


n,Q 4 +n. 


02 04 


0, 


+ n. 


df; 


04 


02 


A 2 = = n , ~ n z O' - !) -Z- 


dQ 


0i 

04 


01 


A 43 = -^ = n y ^ -n z (r- 1) % 

V4 


A 4 


503 

5 /; 


0, 


04 


0+n z (2-y) 
504 01 


45-|^ = ^0 / - 1 ) 

505 


A,, = 


03 04 
01 


0 2 

+ n,— + n,p 

‘ 0i * 


01 ' 5/7 


A> = 


M 

50 5a 


i V 


-«,05 + « v -^- 5 
i41 


0205 , „ 0305 


04 05 


+ n„ _ +n„ _ +n x p—+ n p—+ n.p 


y 0i 


0. 


02 

0, 


.0 

0i 


fit 

0i y 


= 9 


-y 


(Ql +03 2 + 04) 


A„=W =nx 

52 aa 2 

4 d fs 

A,, = — - = n 

3 503 > 

4 d fs 
A SA = ^- = n. 


As = 


504 

M 

505 


05 + Zzl 
01 0 

&±P_S( r _i)& 

a a 

a±Z-f* r -i)& 

a a 

&±il-0( r -i)& 

a a 


yO-n t 


dp _ p y -l 


where = — 

5/7 0, 0! 

specific heat assumption: 


(- 0 5 + 0 2 + 03 + 0 2 ) is derived as shown in Appendix B. For constant 


dp _ y — 1 
dp 2 a 


(e 2 2 +fi 3 2 +a, 2 )=A(' 
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Also, representing the dependent variable vector q as: 


9 = 


Hence, 


P 


<h 


Qx 

u 


<?2 


QJQx 

V 

= 


= 

QjQx 

w 




QJQx 

p 


9s 


9s 


dq 

d Q 


dq 

dQ 


dq x 

dq. 

dq x 

dq x 

dq x 

dQx 

dQ 2 

dQi 

dQ A 

dQ 5 

dq 2 

dq 2 

dq 2 

dq 2 

dq 2 
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dQi 
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dQ A 

dQi 

dq 3 

dq 3 

dq 3 

dq 3 

dq 2 

dQi 

dQi 

dQi 

dQ A 

dQi 

dq 4 

dq A 

dq 4 

dq A 

dq A 

dQi 

dQi 

dQi 

dQ A 

dQi 

dq 5 

dq 5 
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dq 5 
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dQi 
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dp 

dp 


Up 

0 
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0 

Up 
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0 

0 

0 
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0 

0 

0 

0 


— -(/-!)« -(y-l)v -(y-l)w y - 1 


where: 

dq, _ dQ x 
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dq 2 _ dQjQ x 
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da_~ ” ” 


02 

l2 


dQ 2 dQ 4 8Q 5 

-u dq 2 _8Q 2 /Q X _ 


Qi P dQ 2 8Q 2 


Q\ P 


dq 3 _ SQ JQ, _ _ Qi _ - v dq 3 _ dQjQ x 1 _ 1 


dq 2 dq 2 


dq 2 


= 0 


dQi dQ 4 8Q S 
dq 3 _ dq 3 _ 8q 3 _ 


dQi 8Q X 


Qx P dQ 3 8Q 3 Q x p 8Q 2 8Q 4 8Q 5 


0 


NESC Request No.: 07-037-1 



NASA Engineering and Safety Center 
Technical Assessment Report 

Document #: 

NESC-RP- 

07-037 

Version: 

1.0 

Title: 

Reacting Multi-Species Gas Capability for USM3D Flow 

Solver 

Page #: 

52 of 78 


d<? 4 _ dQ 4 /6i _ G 4 _ " vv . Sq 4 _ dQ 4 /Q l 1 _ 1 . dq 4 _ dq 4 _ dq 4 _ Q 

3G, G 1 2 P ’ 3G 4 g'G 4 Gi p’ GQ 2 dQ 3 dQ 5 


And dqs/dQu oqPGQi, GqpGQy, Sq 5 /GQ 4 , Gq 3 /GQ 3 can be obtained as shown in Appendix B and 
shown as the following: 


dq 5 dP y-l 


8Qi op G, 

dq 


p -G 5 +G 2 +G 2 + G 4 2 


y-l 


3G 2 Gi 

^-=-0'- 1 )v=-0'-i)^ 

sg 3 g, 

d( l5 _ i\ t.. i^G 4 


5G_ 

Sg, 


= -(r-l)w = -(/-l)^ 
'di 


5 _ 


5G 5 

Hence, 


dq 

DQ 


y-l 


1 

— u/p 
-vl p 

-wl p 
dp 


0 

Up 

0 

0 


0 

0 

Up 

0 


0 

0 

0 

Up 


— -{y -l)u -(y- l)v — (y — 1) w 

dp 


0 

0 

0 

0 

y-l 


(A.5) 


From Equation 14, the team can obtain: 

dT = -—dp + ^-d p 
p pR 


Hence, 

dT 

dq 


— 0 0 0 ^ 
p pR 


(A.6) 
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Combine Equations (A. 5) and (A. 6), the team can obtain: 


dT 

dQ 


dT dT dT dT dT 


dQ\ dQ 2 dQ, dQ 4 dQ 5 

1 


— 0 0 0 ^ 
p pR 


dT dq 




dq dQ 




0 

0 

0 

0 

Up 

0 

0 

0 

0 

Up 

0 

0 

0 

0 

Up 

0 

(Y~ 1)k 

-O'- l)v 

-(y-l)w 

Y ~ 1 


-v/ p 

— w/ p 
dp 

dp 

~T + 7refdp “/ref (/-!)“ ~YtAY~ l ) V ~ YtAY YtAY 


(A.7) 


p pR dp 


pR 


pR 


pR 


pR 


-T | Y re f (r-ty 2 -YtAY-Vu -YtAY~1)v -Y re[ (Y-l)w YtAY~1) 


p pR 


pR 


pR 


pR 


pR 


The spatial derivatives a dependent variable <j) at the interface between two adjacent cells can be 
calculated as: 

— = r v A^; r= np 

dx n y Ax + n v Av + n_Az 

X V •' z 


*’ = fy =T ’*'’ r ' = 


n„ 


n x Ax + n y Ay + n_Az 
n. 


(/), = — =r_ Acj ) ; r. = — 

dz " ~ n x Ax + n y Ay + n_Az 

and A(j) — </> R — (j) L where (f> R and (/> are the values of </> at the center of right and left. Hence, the 
spatial derivatives of a dependent variable at the cell face can be linearized by the conservative 
Q as: 

r 


d k. =T d(A0 _ r 


dQ 
Therefore, 


dQ 


dQ 


(A.8) 


J 


SAk 

SQr 


=r 


r Wr ^ 
\ d Q K j 
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sq l 


= -r 




Equation (13) can be differentiated with respect to (A.l) to get the viscous flux Jacobian. Each 

dfj 

component of the viscous flux vector can be represented as B jk = and this will result in: 


5 = 


B l2 B 13 B l4 B l5 

B 


B 21 B 22 B 23 B 24 


25 


B 3 \ B 32 B 33 B 34 B 35 

B 41 B 42 B A3 B 44 B 45 

B$ ] B 52 B 53 B 54 B 55 


(A.9) 


Substituting Equations (A. 8), (A. 10), and (A.l 1), each components of the viscous flux Jacobian 
can be expressed as: 

^11 ~ B l2 —B l3 =B [4 = B l5 =0 




Sr. 


Sr,. 


+n„ 


- +n . 


Sr. 


df 2 

8{Qc)j [ 1 x 8(q c )j ,,v s(a) ; '" z s(a),J 


r d K 

^d(Q c ) i 


B *i = 


df; 


Sr,. 


Sr„ 


Sr. 


-+n 


+n_ 


d{Q c )j y d{Q c ) J z d{Q c ) j 


Sr. 


Sr, 


+n. 


-+n_ 


Sr. 


d{Q c )j [ x d{Q c )j y d(Q c ) J z d{Q c )j J 


d/; d(u f /;) s(v/ 3 v ) d(w/;) 
5i d{Q c ) i c(Q c ). S(e c ). S(0 C ). 


+ /c 


sr. 


dT, 


-+n 


+«. 


sr. 


Sfoc) ; y 8{Qc)j z d{Q c )j J 


= u f B 2 j+v f B 3 j + w f B 4j +K 


dT dT dT 

n ,. +n „ — — +/t . 


S((2 C ) ; . y d(Q c ) j z d(Q c )j 


where C represents either cell on the left or cell on the right side of the face and j represents the 
column number. 
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Appendix B. Linearization of Pressure for Real Gases 

For thermally perfect, single-species flow: 

pRT 

P = ~ 

y ref 

Hence, 

dp = ^-dp + ^-dT (B.l) 

7ref y ref 


Since the specific internal energy ( e ) is defined as: 

de= ( ^dT=>—dT= — 
y* f y* f c v 

Also, the specific total energy (E) is defined as: 

E — e + 4 {u 2 + v 2 + w 2 ) 

Hence, 

dE — de+udu + vdv + wdw , or de = dE — (udu + vdv + wdw) 
For a general flow variable ® 

d{p®) = pd® + ®dp => = 

P 


(B.2) 


(B.3) 


Apply the above equation for flow variables, E, u, v, and w to Eq. (B.3), the team can obtain: 

d(pE)- Edp + (w 2 +v 2 + w 2 )dp ~[ud(pu) + vd(pv)+ wd(pw)] 

P 


de = 


(B.4) 


Substitute Equation (B.4) into Equation (B.2): 

1 d(pE)-Edp + (u 2 + v 2 + w 2 )dp-\ud(pu)+ \>d(p\>)+ wd(pw)\ 

y<e{ pc v 

Substitute Equation (B.5) into Equation (B.l), the team can obtain: 
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dp = RTdp-\ [d(pE)-Edp + (u 2 + v 2 + w 2 )dp-[ud{pu) + vd(pv) + wd{pw)§ 


R_ 

c” 


f CJ p 2 , 2 , 2 A 
— E + u +v + w 




dp-ud {pu) -vd(pv) - wd {pw) + d ( pE ) 


y 


= (r-i 


-E + u^ + v 2 + w 


dp-{y -X)ud{pu)-{y -X)vd {pv) 


(r~ ] )p 

-{y-l)wd {pw) + {y-l)d (pE) 

^-dQ, + d Ql + P^clQ, + ^-d Qt + ^JQ, 


dQ x 


dQ 2 


dQ 3 


dQ 4 


32 s 


(B.6) 


where 


Q = 


~q; 


P 

Qi 


pu 

Qi 

= 

pv 

2 4 


pw 

_2 5 . 


_PE_ 


and 


dp 

dQ 

dp 


=(r- 1 ) 


P 77 . 2 . 2 . 2 

— E + u + v +w 




= _p_ + r-i 


y 


a e. 


(- 2 5 + Ql + 2 3 + 04 )- -J- 


dp 


dn =-0'- 1 ) M = -0'- 1 )^- 
d£ 2 2, 


5;? 


~ ■=-0'- 1 ) v =-&'- 1 )% 

d2 3 2i 


dp 


= -(r- 1 V=-(r- 1 )^ 

^2 4 2i 


5/2 


= r-i 


(B.7) 
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Appendix C. Numerical Treatment of Characteristic Boundary 
Conditions in Subroutine UBC2 


Current formulation (constant y = y,)\ 



Figure C-l. Orientation of Boundary Face for Characteristic Variable based Boundary Condition 


Based on the orientation of a boundary face and the nomenclature used in Figure C-l, the 
Riemann invariants can be written as: 


9C =U D +^ = U,+ 2a '' 


^=U L + 


r - 1 

2 a 


L =u b - 


y - 1 

2 a. 


y -1 y — 1 

where U and a are the contravariant velocity and speed of sound, and 


U R = U c + VU ■ As 


u L = u g =u« 


Subscript co denote the freestream condition. Hence, the condition at the boundary can be 
obtained as: 


C„=l(9r+9T); fli =^(9C-9r) 


I *U b +U m > 0 where U m is the contravariant velocity of the moving boundary, then the 
information is traveling from the right-hand-side. Hence, 
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u b= u R +n X U b ~ U r) 
v h =v R +n(U b -U R ) 


+ n z (U b -U R ) 

To obtain the thermal properties at the boundary, isentropic process between R and b is 
assumed, and thus 


pV 


Pr 1 


Pb 


T r/Pr 


r - 1 


and p b = — p b T b 

r 


where 


T b = a 


b , t r = y—\ Pr = p c +vp- a? ; Pr=Pc + vp- as 

Pr 


if u b +u m <0, then the information is traveling from the left-hand-side. Hence, 

U b= U L+ n A U b- U L ) ’ U L= U oo 

v b = v , +n y{U b -U l ) ; v L =v x 
w h = w L+n z (U h -U L ) ; w L =w x 


To obtain the thermal properties at the boundary, isentropic process between L and b is 
assumed, and thus 


Pb 

where 


r - 1 


Pl 1 


Pb = 


TlIP'l 


y - 1 


and p b = — p b T b 

r 


T L = r — ; p l = Poo ; a. = p, 

Pl 


Modified implementation: 

w •=£/„ + ^- = c/ t 

A" 1 

n~=u L +^- = u t 

y L i 


t 2 a 

A -1 
2 A 
A -1 


and 

U r =U c +VU-As; U L =U g =U x ; 


=(«J* =(aL ; A =^o 


(C.l) 


(C.2) 
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where a, is the mass fraction of i-th species and R is the mixture gas constant. Subscript oo 
denotes the freestream condition. Hence, the condition at the boundary can be obtained as: 




(C.3) 


If U b +u m > 0 where U m is the contravariant velocity of the moving boundary, then the 
information is traveling from the right-hand-side. Hence, 

u b = u R +n x (U b -U R ) 

V b = v R+ n y( U b- U R ) 

w b =w R +n z (U b -U R ) 

(a, X = («, )r = ( a s l ’ R b = R R=f((a s ) R ) 


(C.4) 


Once the speed of sound and the mixture gas constant at the boundary are obtained, an 
iterative procedure can be used to determine the temperature and the ratio of specific heats at 
the boundary based on the following correlations. 


a . 


f C ^ 


rM . _ ( c ,l = ( c J R l 

(c p /R {- 1 


i«/ 




f sio \ 

D 

V « Jb 


Yb 


NS 

S = 1 


(c;l 


R 


ajM . 


= 


(C.5) 


(C.6) 


5^=1 

(c^) y/? i( = (a,r + ci^T + + ci^T + u 5 T -\-q,^T -\-ci-iT 


(C.7) 


With the isentropic assumption between the upwind cell and the boundary, the mixture 
density and pressure at the boundary can be calculated based on the temperature (77) and 
mixture density ( p c ) of the upwind cell as: 


P b = 




1 

n 1 1 


and p b = —p b R b T t 

Yrel 


(C.8) 


Total energy at the boundary is then calculated from temperature, species mass fraction, and 
velocity components at the boundary based on the caloric equation of state. 
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If U b +u m <0, then the information is traveling from the left-hand-side. Hence, 

u b =u L +n x (U b -U L ) ; u, = u x 
v b =v L +n y (U b -U L ) ; v, = v a 
w b = w L + n z (U b -U L ) ; w L =w x 
(«, )/, = («, )l = («, L ; R h = R R = k 

The procedure stated above used to solve Equation (C.5) through (C.7) can be used to obtain 
the thermal properties at the boundary. Finally, the total enthalpy at the boundary can be 
calculated as: 

e „=£ (c p - R \ dT + ik+v 2 b +w;) and E g = E x (C.10) 
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Appendix D. Numerical Treatment of Viscous Wall Boundary 

Conditions 



Figure D-l. Orientation of a Wall Boundary Face 


Based on the notations used in Figure D-l, the boundary conditions for a solid wall boundary can 
be summarized as: 


u b = «, ; v b = v w ; W b =w w 

u g = 2 u b -u c \ v g = 2v b - v c ; = 2 w b - w c 

dp 

— = 0 => p b = p c +Vp-As ; p = 2p b - p c 
on 


Current formulation: 


1) adiabatic wall: 


^ = 0 
dn 


T b =T c - T=T C 


Pb = Xref 




2) isothermal wall: 

Tb = p h = /jef — - 

‘ b 

P g =2p b -Pc 
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Modified implementation: 

dec 

-^ = o^ («,)*=(«,)*; («,),= («Jc 
on 

1) adiabatic wall: 


T 

— = 0 


T b = T c ; r g = r e 
Pb . _ JT 

6 ref ’ e b Q Q dT 

1 b 


— £. * : 

g ref rri ’ 

1 g 


C, dT 


and E h = e b + i(u 2 + v 2 + w 2 ) ; £ g = e g + \ (w 2 + v 2 + w 2 ) 


2) isothermal wall: 


T,=T => 


E b = \o C v dT+ 2( U l + V b+ W b) and A,= 


^ ref 


Pb 


P g = 2 Pb-Pc 


T,=kA and E,=fc,dr + i(« g 2 +^+^) 
Pg 
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Appendix E. Description of USM3D Input File 


USM3D (v6.0) - Sample Input File 


Mach 

alpha 

beta 

ReUe,mil 

Tinf , dR 

itwall 

Tw/Tinf 

ipwall 

0.84000 

3.06000 

0.00000 

11 . 7 

460.00000 

0 

-1.0 

0 

sref 

cref 

bref 

xmc 

ymc 

zmc 



0.52550 

0.66700 

1.00000 

0.00000 

0.00000 

0.00000 



ioverset 

impl 

cf 11 

iramp 

cf 12 

cf lmin 

GS tol 

crelax 

0 

1 

-150.000 

50 

150. 

5.0 

-20 

1.00 

itimeacc 

delta t 

ntstep 

res step 

imvgrd 

isolavg 

timebgn 


-2 

0.3 

500 

-3.0 

1 

0 

1500.0 


irest 

mstage 

iresmth 

dqmax 

p break 

pmin 

limiter 


0 

3 

1 

0.5 

0.050 

0.001 

0 . 


nupdate 

nwrest 

nwf lo 

nwf lobgn 

ipltqn 

idiagnos 

nodeypl 


50 

75 

100 

1500 

2 

0 

0 


iorder 

lapl-avg 

high-bc 

ifds 

ivisc 

itr ip 

ev lim 


2 

1 

1 

1 

3 

0 

0.0 


ncyc 

nengines 

nsinkbc 

nrotor 

compF&M 

p bcl002 

cldes 


10 

1 

1 

1 

0 

0.0 

0.0 


jetl 

fuel 

gamma j 

pjet 

pO j et 

Rratio 

TO j et 


103 

0.02000 

1.40000 

0 .71430 

1 .42857 

0.77400 

3.17390 


jetl 

vec (1) 

vec (2) 

vec (3) 

engsim 

nstation 

idirec 


103 

1.00000 

0.00000 

0.00000 

0 

10 

1 


jetl 

fuel 

gamma j 

pjet 

pO j et 

Rratio 

TO j et 


102 

0.02000 

1.40000 

0 .71430 

1 .42857 

0.77400 

3.17390 


jetl 

vec (1) 

vec (2) 

vec (3) 

engsim 

nstation 

idirec 


102 

1.00000 

0.00000 

0.00000 

1 

10 

1 


jetl 

inner R 

outer R 

cen X 

cen Y 

cen Z 



102 

1 . 0114 

1.39846 

23.646 

8.713 

-3.322 



(r-inner R) / (outer R- 

-inner R) 






0.168,0.350, 

0. 500, 0. i 

600,0.700, 0 

. 800, 0.900 

,0.950,0.975,1.00 



pO jet 








1.07143 1.07143 1.07143 1.07143 

1.07143 1 

.07143 1.07143 1.07143 1.07143 

1.0714. 

TO jet 








1.32610 1.32610 1.32610 1.32610 

1.32610 1 

.32610 1.32610 1.32610 1.32610 

1.326H 

V circular /Vmag 







0.0168 0.035 

0.050 0 

.060 0.070 

0.080 0.090 0.095 0 

.0975 0.1 



V radial /Vmag 







O 

o 

o 

o 

0 . 0 . 0 . 

0 . 0 . 0 . 






sinkid 

outer R 

cen-X 

cen-Y 

cen-Z 

massf lx 

omega 


111 

1.0 

0.0 

0.0 

0.0 

3.0 

-0.5 


rotid 

inner R 

outer R 

cen-X 

cen-Y 

cen-Z 



502 

0.00000 

24 . 0000 

0.0 

0.0 

0.0 



Advratio 

thrst c 

torq c 

irotat 

rotsim 




0.20 

0.00861 

0.00 

1 

0 




Xrcg 

Yrcg 

Zrcg 

ref len 

iroll 

ipitch 

iyaw 


0.0 

0.0 

0.0 

1.0 

0 

1 

0 


ang mean 

ang ampl 

red freq 

idef ormg 

imvMOMCEN 




0.00 

2 . 41 

0.0808 

0 

0 




ikeord 

icons 

nstagek 

t dtfact 

t intsity 

mut/mul 

ratiokp 

dkemax 

1 

1 

5 

1.0 

l.e-3 

1.0 

0.00 

0.25 

ini 

ilhg 

iwallf 

icompCorr 

itempCorr 

itk 

isk idt proc 

7 

-14 

0 

0 

0 

2 

1 

0 
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flkemax <=== THIS IS A TEMPORARY VARIABLE FOR CAPPING MAX fl in LB MODEL 

1.0 

$SPECIES 

H20 02 CO C02 C7H8 H2 0 H OH N2 

$REACTIONS 
13 


REACTION :H20, 02, CO, C02 , 

C7H8 

CM 

X 

0, 

H, 

OH, 

N2 

1, 

4.4963E9, 1.00, 2.679E4 

, o, 

2, 






0 . , -3.5, 7 . , 0 . , 

-1 . 

, 4., 

0. , 

0., 

0. , 

0. 


0. , -1.0, 0. , 0. , 

-0.5 

, 0., 

0. , 

0., 

0. , 

0. 

2, 

1 . 7000E13, 0.00,24070., 

0, 

0, 






0., -1., 0., 0., 

0. 

, -1-, 

0. , 

0., 

2., 

0. 

3, 

2 . 1 900E13 , 0.00, 2590., 

0, 

0, 






o 

o 

o 

V — 1 

0. 

, -1., 

0. , 

1., 

-1. , 

0. 

4, 

6 . 0230E12 , 0.00, 550., 

0, 

0, 






O 

o 

o 
' — 1 

0. 

, 0., 

1. , 

0., 

-2., 

0. 

5, 

1 . 8000E10, 1.00, 4480., 

0, 

0, 






o 

o 

o 

o 

0. 

, -1., 

-1. , 

1., 

1. , 

0. 

6, 

1 . 2200E17 , -0.91, 8369., 

0, 

0, 






o 

o 

\ — 1 

1 

o 

0. 

, 0., 

1. , 

-1., 

1-, 

0. 

7, 

4 . 0000E12 , 0.00,10030., 

0, 

0, 






0., 0., -1., 1., 

0. 

, 0., 

0. , 

1., 

-1. , 

0. 

8, 

3 . 0000E12 , 0.00,2. 50E4 , 

0, 

0, 






0., -1., -1., 1., 

0. 

, 0., 

1. , 

0., 

0., 

0. 

9, 

1 . 0000E16, 0.00, 0., 

999, 

0, 






o 

o 

o 

o 

0. 

, 0., 

-1. , 

-1., 

1. , 

0. 

o 
\ — 1 

2 . 5500E18 , -1.00,59390., 

999, 

0, 






0., 1., 0., 0., 

0. 

, 0., 

-2. , 

0., 

0., 

0. 

\ — 1 
\ — 1 

5 . 0000E15, 0.00, 0., 

999, 

0, 






o 

o 

o 

o 

0. 

, 1., 

0. , 

-2., 

0. , 

0. 

CM 
\ — 1 

8 . 4000E21 , -2.00, 0., 

999, 

0, 






o 

o 

o 
' — 1 

0 . 

, 0 ., 

0 . , 

-1., 

-1., 

0 . 

13, 

6 . 0000E13, 0.00, 0., 

999, 

0 , 






0 ., 0 ., -1., 1., 

0 . 

, 0 ., 

-1. , 

0 ., 

0 . , 

0 . 


c=========== 

C LINE 1 
c Title 

C ===========' 

c LINE 2 
c xmach 

c 
c 
c 

c alpha 

c beta 

c reue 

c tint 

c itwall 

c 
c 

c twtinf 

c 

c 


- free form upto 80 characters 


- freestream Mach number 

for xmach < 0, velocity will be initialized by 

a fixed mach value of 0.4. Some usefulness for supersonic cases, 
however *not* recommended. 

- angle of attack, deg. (in X-Z plane) 

- side-slip angle, deg. (in X-Y plane) 

- freestream Reynolds number per unit length (millions) 

- freestream temperature (degrees Rankine) 

- wall temperature boundary condition flag 
= 0 adiabatic wall temperature 

= 1 specified wall temperature 

- wall temperature (temperature at wall divided by temperature 
of freestream) 

if twtinf<=0, twtinf taken as freestream stagnation temperature) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ipwall - type of extrapolation of pressure to wall 
= 0 zeroth order 

= 1 first order (linear extrapolation) 


LINE 2 

igmodel - flag to specify gas model 
= 0 Air (calorically perfect) 

= 1 Any gas with specified Sutherland const and gamma 
(calorically perfect) 

= 2,-2 Any gas with specified Sutherland const, gamma, and 
gas name (to set temperature variation of gamma) 

= 3,-3 Two-gas model with specified gamma, Sutherland const, 
stagnation enthalpies, and gas names 
-2 or -3 is Recommended 
gammal , gamma2 - specific heat ratios (gamma) 
htotl,htot2 - total enthalpies (non-dimensional) 

numspec - total number of species constituting 1- or 2-gas system 

NOTE : For igmodel values of 2 or 3, user must have "thermo. inp" 

file of McBride and Gordon. 

: For igmodel=2 or 3, you can set gammal<=0; code will 
calculate the value internally based on freestream temp. 

: For igmodel=3, you can set Htotl<= 0; code will calculate 
the value internally. 

: For igmodel=3, and exit BC -102, you can set Htot2<= 0; 
code will calculate the value internally. 


c LINE 3 


c gasname - Name of a gas (must be part of McBride-Gordon table) 

c ispec_grp - specie group index, must be 1 or 2 (2-gas model) 

c mole_frac - Mole fraction of the gas in a given specie group 

c 

c LINE 5 

c sref - reference area (half-area for half configuration) 

c cref - reference length 

c bref - reference span 

c xmc - moment center in x-direction 

c ymc - moment center in y-direction 

c zmc - moment center in z-direction 

c 

c 

c LINE 6 
c ioverset 
c 
c 

c impl 

c 
c 
c 
c 
c 
c 
c 
c 


- overset grid flag 
= 0 single grid 

= 1 overset grid 

- implicit time stepping flag 

= 0 explicit Runge-Kutta time stepping 
= 1 implicit Gauss-Seidel time stepping 

= 1 UPDATQ, timestep scaling; UPDATQ enforces positive- 

definite updates for density & internal energy 
= 10 UPDATQ, timestep scaling off 

= 2 UPDATQ1, timestep scaling; UPDATQ1 *does not* 

change solver-provided solution updates 
= 20 UPDATQ1, timestep scaling off 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C LINE 7 

c itimeacc 

c 

c 

c 

c 

c 

c 

c delta_t 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c ntstep 

c 

c res_step 

c 

c imvgrd 

c 

c isolavg 

c 
c 
c 


cf 11 

iramp 
cf 12 
cf lmin 
GS tol 


crelax 


= 3 UPDATQ2, timestep scaling 

= 30 UPDATQ2, timestep scaling off 

code will internally reset impl to 1 after raising 

appropriate flags. 

time step or initial elf number 

< 0 local time stepping or ramping 

number of cycles for ramping cfl number 

final value of cfl number for ramping 

minimum allowable CFL for variable time step strategy 
Tolerance criterion for automatic termination of Gauss-Seidel 
iteration process termination for Flow and S-A equations 
recommended value ( l-dq_prev/dq_new) is 0.001-0.005. 

Negative value will be interpreted as G-S stages 
by the code for solving linear problem for mean flow. 

Minimum threshold value for stages is set to 5. 

Under /over-relaxation coefficient for Newton's scheme. 

0.7 as recommended by Dr. Paul Pao. 

Without this relaxation, Newton's scheme produced jagged 
hysteresis of Lift and moment coefficients for 
NACA 0012 viscous pitching case. 

*CAUTION* itimeacc=0, 1 or 2 options also will use relaxation 
coefficient. For itimeacc=2, relaxation is not typically 
needed. Therefore it is important to set it to 1.0 when 
not using Newton's scheme ( itimeacc=-2 ) . 


0 Not time accurate; assumes steady state final solution 

1 Time-accurate solution with subiterations using 
Crank-Nicholson 

2 Time-accurate solution with subiterations using 
3-point backward differencing and pseudo-time 

-2 Time-accurate solution with subiterations using 

3-point backward differencing without pseudo-time (Newton) 
time step increment (nondimensional ) 

HOW DOES ONE DETERMINE VALUE FOR TIMESTEP? 

-> Think of it as how many timesteps one needs for a fluid 
particle to travel some characteristic length, Lchar, 
e.g. wing chord, at freestream velocity. A rule of thumb 
from Jim Forsythe, Cobalt Solutions, is N=50 steps. 

Scott Morton of the Air Force Academy suggests between 
100 and 1000 steps. 

-> For USM3D nondimensionalization, the relationship for 
timestep and number of steps, N, to traverse Lchar: 
delta_t = Lchar/xmach/N 
number of timestep increments 

(Note that total "iterations" will be approx. ntstep*ncyc) 

level of subiteration convergence before taking next step 

Typical value might be -2 or -3. 

invokes prescribed grid motion 

(additional input required, see LINE 5a below) 

solution averaging flag 

0, turn-off solution averaging (for steady-state, it will 

be turned off by the code) 

1, perform time averaging for mean flow solution qunatities 


NESC Request No.: 07-037-1 



NASA Engineering and Safety Center 
Technical Assessment Report 

Document #: 

NESC-RP- 

07-037 

Version: 

1.0 

Title: 

Reacting Multi-Species Gas Capability for USM3D Flow 

Solver 

Page #: 

67 of 78 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


(extra memory equivalent to 5*nnode) . Also, mean values 
of aerodynamic quantities will be printed. 

= 2, In addition to 1 above, perform solution averaging for 
turbulent-stress-related qunatities (extra memory 
equivalent to ll*nodes) 

timebgn - beginning timestep for initiating time-averaging of solution 
(based on the "ntstep" timestep numbering, not "ncyc") 

NOTES: 1) "ncyc" becomes the number of subiterations+1 used to converge 
solution to Q_(n+1) 

2) Need to monitor hist.subit for subiteration convergence 


c LINE 8 


c irest = 

c = 

c mstage 

c 
c 

c iresmth - 

c = 

c = 

c 

c 

c 

c 

c 

c 

c 

c 

c dqmax 

c 

c = 

c = 

c pb r k = 

c pminimum= 

c 
c 

c bcoef 

c = 

c = 

c = 

c = 

c 

c 


0 no restart 

1 restart 

number of Runge-Kutta stages (for explicit scheme) 

For implicit scheme mstage < 0 specifies interval at which 
all the cells are reset as second-order cells 
implicit residual smoothing (for explicit scheme) 

1 smoothing on every m-stage cycle 
-1 smoothing only on odd m-stage cycles 
For implicit scheme iresmth specifies technique to mark 
first order cells 

= 0 => do not set cells to first order locally 
= 1 => use p<=pminimum or rho<=pminimum switch to set 
first order cells 

= 2 => use p<=pminimum or rho<=pminimum or pt>=1.2*pt0 
switch to set first order cells iresmth=2 will 
only be used for non-engine, non-rotor cases 
maximum allowable (-dp/p, -drho/rho) in solution 
(used to scale CFL number) 

0.5 suggested 

greater than zero --> No Limiter Applied 

threshold pressure for lowering local CFL toward CFLMIN 

Minimum allowable level of non-dimensional pressure and 

threshold for setting cells to 1st order 

<0 also resets Ist-order flagged cells to 2nd-order 

Limiter coefficent 

0. no limiting 

1 . MINMOD limiter 

2 . SUPERBEE limiter with variable compression factor 
-2. will involke automatic selection of limiter 


c LINE 9 


c 

c 

c 

c 

c 

c 

c 

c 

c 


nupdate 

nwrest 


nwf lo 


number of iterations between updates of the time step 
(not used for time accurate computations) 
number of iterations between updates of the binary 
restart file 

for UNSTEADY runs ( itimeacc . NE . 0 ) , this is restart file 
updates are based on the number of "ntstep", not "ncyc" 
number of iterations between writing of instantaneous 
NEW flow-file (pr j_name . f lo . iter # ) . Typically used for 
time-accurate runs . 


NESC Request No.: 07-037-1 



NASA Engineering and Safety Center 
Technical Assessment Report 

Document #: 

NESC-RP- 

07-037 

Version: 

1.0 

Title: 

Reacting Multi-Species Gas Capability for USM3D Flow 

Solver 

Page #: 

68 of 78 


nwflobgn - number of iterations beyond which writing of instantaneous 
NEW flow-file (pr j_name . f lo . iter # ) will be triggered. 

The parameters nwflo and nwflobgn relate to outer loop 
time steps (ntstep) for an unsteady flow analysis 
and to inner loop time steps (ncyc) for steady flow analysis, 
ipltqn = 0 suppress output of post analysis of flow field file 

= 1 write unformatted q-node file for post analysis of flow field 

= 2 write formatted q-node file for post analysis of flow field 

= 3 write formatted q-node file for postprocessing by TECPLOT 
= 4 write binary TECPLOT file (NOT CODED YET) 

= 5 write unformatted q-node file for postprocessing by FIELDVIEW 

X =13 write unformatted q-node file for using Tecplot converter 

X =23 write formatted q-node file for using Tecplot converter 

=15 write unformatted q-node file for using FIELDVIEW converter 
=25 write formatted q-node file for using FIELDVIEW converter 
=-25 write solution restart file for using FIELDVIEW converter 
also write formatted q-node file and do not delete 
X =16 write unformatted q-node file for using Ensight converter 

X =26 write formatted q-node file for using Ensight converter 

X = **indicates that this option has not been tested and** 

**not recommended for use.** 
idiagnos= 0 no output 

1 write history of minimum static pressure, maximum 
Mach number, maximum eddy viscosity in info.diag file. 

2 write a flow file with 5 auxiliary variables in place of Qs 
for unsteady flow simulations. These variables are, 
instantaneous and time-average values of vorticity, 

eddy viscosity, and local CFL number. 

3 write info.diag2 file in Tecplot format that contains 
several variables at cell-centers for the cells associated 
with user-specified boundary faces. These boundary faces 
must be supplied in project. trip file. You will need to 
create a project. trip file for idiagnos != 3 as well, if 
itrip != 0 (below) . NOTE: Diagnostics related to 
idiagnos = 3 is available only for ivisc > 0. 

4 write info. probe. N file in Tecplot format that contains 
history of Qs for each probe specified in pro j ect . probe 
file. N will be substituted by the probe number. 

5 write "Sol .Volume . USM3D" file in PLOT3D format that contains 
Qs at desired locations in the grid as read from the file 
"Grid. Volume . TEACC" in a pre-processing step. In addition, 
"Sol . Exit . TEACC" file will be read at appropriate time to 
derive exit pressures for appropriate locations on the 
boundary. These locations are written out in the file 
"Grid. Exit. USM3D" . 

nodeypl= # Boundary node number where output of u+/y+ is desired 
0 no output 


c LINE 10 

c iorder 

c 

c 

c 

c 


= order of flux computation 

= <0 (INTEGER only) Automatic selection of order based on 

prescribed magnitude of mean-flow residual error drop 
= 0 Automatic selection of order based on residual error drop 

by one order 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


laplavg 


ihibc 

ifds 


ivisc 


itrip 


ev lim 


1 first order 

2 Taylor series with full gradient computation 

0 Inverse-distance averaging to all nodes (Vers. 1.0) 

1 Pseudo-Laplacian averaging to all nodes 

-1 (Ps. Lapl . on interior, Inv. Dist. on boundary nodes) 

0 First-order accurate boundary conditions (Vers. 1.0) 

1 Second-order accurate boundary conditions 
spatial differencing parameter for Euler fluxes 

0 flux-vector splitting 

1 flux-difference splitting (Roe scheme) (recommended) 

2 AUSM+ (need to evaluate for USM3D (as of 06/17/04). 

3 HLLC (need to evaluate, only for stationary cases 

as of 07/12/05) 

viscous/inviscid interaction flag 

0 inviscid 

1 laminar (not active in version 2.0) 

2 Spalart-Almaras, full viscous 

-2 Spalart-Almaras, full viscous with DES 

3 Spalart-Almaras with wall function 

-3 Spalart-Almaras with wall function with DES 

6 K-Epsilon turbulence model 

-6 K-Epsilon turbulence model with DES 

-60 K-Epsilon turbulence model with Multiscale capability 

7 K-Epsilon turbulence model, Carlson's formulation 

8 SST turbulence model 

-8 SST turbulence model with DES capability 

-80 SST turbulence model with Multiscale capability 
flag to activate flow tripping 

0 no tripping 

1 activates tripping (useful mostly for k-e turbulence model 
when used for low Reynolds number flows) 

-1 enforces laminar flow region (useful mostly for SA and 
SST model when a certain flow region has to be forced 
to a laminar flow state. These two models generally 
provide a fully turbulent flow on the body) . 

NOTE: for itrip value other than 0, user will need to 

(1) define trip patches in the existing grid 

using PREDISC program (contact: Richard Campbell, 
LaRC) 

(2) create a file project. trip and list all the 
trip patches and appropriate trip values. 

Refer to the format of project. trip file on 
the USM3D website. You will need to create 
project. trip file for itrip = 0 as well, if 
idiagnos = 3 (above) . 

Eigenvalue limiter value 


NOTE: iorder automatic switch should be used only for steady-state cases. 
If you want to force first order computations based on 
previous second-order solution upon restart, you must 
set iorder = 1 . Negative iorder value will not let you 
restart the solution in first order mode. 


c LINE 11 
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c 

c ncyc 

c 

c ngine 

c 

c nsinkbc 

c 

c 

c nrotor 

c 

c compF&M 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c p_bcl 002 

c 

c 

c cldes 

c 

c 

c 

c LINE 11a 


STEADY-STATE PROBLEMS: number of time step cycles 

UNSTEADY PROBLEMS: number of subiterations to "ntstep" 

number of jet engines (maximum of 4) pmh 

(additional input required, see LINE 9a below) 
number of special intakes with specified mass flux 
required be type in .mapbc: 111,211,311,411. 

(additional input required, see LINE 9b below) 
number of rotors/propellers 

(additional input required, see LINE 9c below) 

0 do not include "comp.FandM" file; no special integration 

1 integrate component F&M's over patches prescribed in 

file "comp.FandM' at the end of a run. write on to tet.out 
-1 integrate component F&M's over patches prescribed in 
file "comp.FandM' continuosly during run. write the 
CL, CD, CDv, CM, and 6 components of body forces and 
moment history on "pro j ect . FMhistCOMPN" file where, 

N stands for component index. You will get as many files 
as the number of components in "comp.FandM" file. 

Also, write the last iteration F&M information on to 
tet.out file 

2 automatically treat every existing patch in the grid as a 
component and write forces and moments on to tet.out 

at the end of a run. "comp.FandM" file not read/needed. 

-2 automatically treat every existing patch in the grid as a 
component and write history of forces and moments in 
the seprate files. You will get as many files as the 
number of active patches in the grid. Also, write 
last iteration F&M information on to tet.out. "comp.FandM" 
file not read/needed. 

3 merger of features activated by compF&M=l and compF&M=2 

-3 merger of features activated by compF&M=-l and compF&M=-2 

-12 merger of features activated by compF&M=-l and compF&M=2 

-13 merger of features activated by compF&M=-l and compF&M=3 

non-dimensional plenum pressure for the BC 1002 

<= 0.0 Use freestream pressure (1/gamma = 0.714) 

> 0.0 Use the specified value 

prescribed CL for matching lift coefficient via feedback loop 
(cldes=0. means no feedback and maintain prescribed alpha) 


-1-k-k-k'k'k-k-k'k'k-k-k'k-k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k-k-k-k'k'k-k-k'k'k-k'k'k 

;*** NOTE: These lines are included ONLY IF ngine. ne.O *** 

-i'k'k'k-k-k'k'k-k-k'k'k'k'k'k'k'k'k'k'k-k'k'k'k-k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k'k'k-k-k'k'k'k'k'k'k-k 


ibetype 


boundary condition flag 

Second jet or fan exhaust: 
Primary jet exhaust by: 


betype 

betype 


n*100+3 (optional) 
n*100+2 


where "n" is the engine number 


c 
c 
c 
c 
c 

cvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 
c RULE: If second jet or fan exhaust is utilized, it's 

c characteristics must be prescribed before those 

c for the primary jet b.c. 

^-.AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
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c fuel 

c gammaj 

c pj et 

c pOjet 

c Rratio 

c TOjet 

c vec(i) 

c engsim 

c 
c 
c 

c nstation 

c 

c idirec 

c 

c 

c 

C NPR 

C TTR 

C RMACHJ 

C 

c 


= fuel fraction - in % of total exhaust mass flow 
= ratio of specific heats for (hot) jet 
= static nozzle pressure (usually set to p_infinity) 

= stagnation pressure of jet 

= spec, gas constant of air/spec. gas constant of jet 
= stagnation temperature of jet 

= directional cosines of exhaust jet (i=l,2,3 for x, y, z) 

= 0 uniform exhaust properties such as pOjet, tOjet 
= 1 radially varying exhaust properties 

will need radial distribution of total pressure, 
total tempreature, tangential and radial velocities. 

= number of radial locations at which exhaust properties 

specified (only useful for engsim=l) . LIMIT: nstation <=20. 

= direction of tangential velocity (only useful for engsim=l) 

= 1 clockwise swirl on engine exhaust looking downstream 

= -1 counter-clockwise swirl on engine exhaust looking downstream 

= nozzle pressure ratio, pOjet/pjet 
= total temperature ratio, TOjet/TOinf 
= exit mach number of jet 


LINE lib 

sinkid 


outer^R 
cen-X 
cen-Y 
cen-Z 
ma s s f 1 x 


omega 


sink be identification 

sink 1 (BC ill), sink 2 (BC 211), sink 3 (BC 311)... 

#of sink bes specified earlier in the input as nsinkbc, 

nsinkbc <= 4 

outer radius of sink 


sink BC ' s center coordinates 

mass flow rate through the sink (unit ???) 

-> nondimensional_massflx = dimensional {mdot/ ( rhoinf*ainf) } 
ratio of fan-tip velocity to freestream velocity 
positive for clockwise looking downstream 
negative for counter-clockwise looking downstream 


c 

LINE lie 



c 

rotid = rotor 

identification 

c 

rotor 

1 - BC 

501,502 

c 

rotor 

2 - BC 

601,602 

c 

rotor 

3 - BC 

701,702 

c 

rotor 

4 - BC 

801,802 

c 

inner R = inner 

radius 

of the rotor 

c 

outer R = outer 

radius 

of the rotor 

c 

cen-X, cen-Y, cen-Z = 

rotor i 

center coordinates 


Advratio 

thrst_c 

torq_c 

irotat 


rotsim 


= advance ratio (free stream velocity/rotor tip speed) 

= thrust coefficient of rotor 
= torque coefficient of rotor 

= direction of rotation (viewing from front/top) 

= 1 clockwise rotation 
= -1 counterclockwise rotation 

= 0 rotor be with uniform jump in pressure and/or swirl 
= 1 rotor be with pressure and/or swirl jump as f(r) 
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c 

c LINE 12a 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


.CAUTION: These two lines of input read/needed ONLY IF imvgrd != 0 in LINE 4 
.CAUTION: Current implementation set up to allow only 1 DOF for rotation. 


xrcg, yrcg, zrcg 
ref len 

iroll 

ipitch 

iyaw 


angmean 
angampl 
rf req 


idef ormg 


imvMOMCEN 


origin of the rotation axis 

reference length, used for calculating reduced 
frequency, based on grid coordinates 

grid rotation about x-axis (rotation axis | | x-axis) 

1 for activating rolling motion 

0 for NOT activating rolling motion 

grid rotation about y-axis (rotation axis | | y-axis) 

1 for activating pitching motion 

0 for NOT activating pitching motion 

grid rotation about z-axis (rotation axis | | z-axis) 

1 for activating yawing motion 

0 for NOT activating yawing motion 
mean angle of rotation (in degrees) 
amplitude of rotation (in degrees) 

reduced frequency of oscillations (rfreq=k/Lchar) where 
k=2 *pi* f *Lchar/ (2*V_inf) , where f is frequency, Lchar is 
characteristic length, where all quantities are in 
dimensional form 

The number of "steps" required to cycle through one 
2*pi period of motion can be computed via. 

Nperiod = pi*Lchar/ ( k*xmach*delta_t ) 

grid deformation flag, important only for itimeacc=2. 

0 for solid-body rotation of grid, 

value of 0 will switch off Geometric Conservation Law. 

1 for grid deformation to accommodate body motion, 

0 do not rotate moment refernce center with the grid 

1 rotate moment refernce center with the grid 
( recommended) 


cNOTE 

c.. angle for pitch (positive angle for nose up) 
c.. angle for roll (positive angle for right wing up) 

c.. angle psi for (positive angle for left turn, wind blowing in pilot's 
c . . left ear) 

c 

NOTE: VGRID generated proj.cogsg file will now be protected. (Aug 31, 2005) 

c 


LINE 13a 
..CAUTION: These 

ikeord 


three lines of input read/needed ONLY IF ivisc = 6, 7, or, 8 

- spatial order of accuracy for 2-equation 
turbulence models 

= 1, first order 
= 2, second order 
Recommended value is 1 

- Flag to activate either incompressible or compressible 
formulations of for 2-equation turbulence models 

= -1, incompressible 

= 1, compressible, conservative formulation 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


nstagek 
t dtfact 


t_intsity 

mut/mul 


ratiokp 


dkemax 


c LINE 13b 

c ini 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c ilhg 

c 

c 

c 

c 

c 

c 

c iwallf 

c 
c 
c 

c icompCorr 

c 

c 


Recommended value is 1 

- Gauss-seidel subiterations for the implicit solution 
of 2-equation turbulence models 

Recommended value is 10 

- parameter to decrease/increse CFL number for the 
solution of 2-equation turbulence models 

> 1.0 decreases CFL number 
< 1.0 increases CFL number 
Recommended value is 1.0 

- freestream turbulence intensity 
Recommended value is 1.0e-3 

- ratio of freestream turbulence viscosity to the 
freestream laminar viscosity 

*NOTE* : You can specify either turbulence intensity 
(t_intsity) < 0 or ratio of turbulence-to-laminar 
viscosity (mut/mul) < 0 to activate USM3D's default 
values of freestream turbulence intensity 
[=sqrt ( 9 . e-9 ) ] and mut/mul [=0 . 00 9 ] 

- Upper limit on turbulence kinetic energy in relation 
to static pressure 

Recommended value: 0.0 for first attempt at solution 
If convergence problems are encountered specify 
a value in the range of 0.05-0.1 

- Maximum allowed update of turbulence model variables 
at any given time step (=Delta (var) /var) 

Recommended value: 0.25 


- Integer to activate different formulations of two 
equation turbulence models 

0 = Basic turbulence stress model, assumes that 

stresses are linearly related to the strain 

1 = Nonlinear model of Shih, Zhu, and Lumley 

Ref. AIAA- 95-224 6 
7 = Nonlinear model of Girimaj i 
Ref. NASA CR-198243, 1995 

*NOTE* : Currently, only k-eps model has provision to 
activate different formulations. For SST model, only 
one value (inl=0) is allowed and user input will be 
overridden if an inconsistent value is specified for 
this model 

- Near-wall modification of k-eps turbulence model 
-14 = Formulation proposed by Launder and Sharma 

-3 = Formulation proposed by Launder and Sharma with 
wall function 

+3 = No near-wall modification of k-eps turbulence 

model but use wall function 
Recommended value: -14 

- Flag to specify wall function grid 

0 = grid suitable for integration all the way to wall 

1 = wall function grid 
Recommended value: 0 

- Formulations for correction to the turbulence models 
related to compressibility at high speeds 

0 = no correction applied 
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c 

c 

c 

c 

c 

c 

c itempCorr 

c 

c 

c 

c 

c 

c itk 

c 

c 

c 

c 

c 

c 

c 

c 

c isk 

c 

c 

c 

c 

c idt_proc 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c LINE 13c 


1 = Sarkar ' s formulation. Ref: J. of Fluid Mechanics, 

Vol. 227, 1991 

2 = Zeman's formulation. Ref: Physics of Fluids, 

Vol. 2, 1990 

5 = Wilcox's formulation. Ref: AIAA Paper 92-0439 
Recommended value: 1 

- Flag to adjust (for linear k-e model only) eddy 
viscosity for high-temperature flowfields, 

for example, hot jets. 

0 = no adjustment made 

1 = temperature correction model of Hamid 

Ref: J. of Fluids Engg . , Vol. 126, No. 4, 2004. 

- Strategies for applying wall boundary condition for 
the turbulence dissipation variable of k-e model 

0 = isotropic component of dissipation set to zero 

(Launder and Sharma model) 

1 = zero gradient of dissipation rate at the wall 

(Lam and Bremhorst model) 

2 = dissipation proportional to turbulence kinetic 

energy grdient at the wall 
Recommended value: 2 

- Flag to apply correction from turbulence kintetic 
energy equation to the mean-flow energy equation 

0 = no modification 

1 = apply corrections 
Recommended value: 0 

- Flag to activate a specific procedure to compute local 
time-steps for steady-state solutions 

0 = use face-based subroutine (UCTIME_MOD) 

1 = use node-based subroutine (UCTIME_KE) 

For SA and SST models, idt_proc is reset to 0 
internally . 

NOTE: This input variable is introduced in light of 
the observed sensitivity of a CEV solution (Run 314) 
to the differences in the procedure used to compute 
local time steps. This option may be eventually 
dropped from the export version 6.0. 


c flkemax - Maximum allowed value for function fl in dissipation 

c rate equation for Lam Bremhorst model (itk=l) 

c Recommended value : 1.0 

c for itk = 0 and 2 this input variable is not used. 


Supplemental Description of Input for Real Gas and Chemical Reactions 


$SPECIES: flag to activate real-gas model for single or multiple chemical species 


The above line is followed by a list of species names involved in this test case. The name of 
each species can be either separated by blank space(s) or comma. An example is shown below: 
H 2 0 0 2 CO C0 2 C7H8 H 2 O H OH N 2 


$REACTIONS: flag to activate chemical reactions 
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The above line is followed by the number of chemical reaction steps (NR). If NR is negative, 
then equilibrium chemistry will be used. If the specified value is greater than 0, then the finite- 
rate chemistry will be used where the value indicates the number of reaction steps. The format 
of finite-rate chemistry model is listed below. As shown in Equation (E.l), a general system of 
chemical reactions can be written in terms of its stoichiometric coefficients: 


( v'ij and Vjj )and the z'-th chemical species ( ® k and O'f ) for the j-th reaction: 


NS Kr ■ NS 

TV® Y/.o 

Lu Jl l Tj r Jl l 


i = 1 


K 


b,j '- 1 


(E.l) 


The source term of z-th species ( co i ), due to chemical reaction, can be evaluated as: 


co. 


A, ii 


f v* 

pa, 


i=i 


M, 


-K, n 


v ‘ i y 


i=\ 


r \ r f 

pa, 


(E.2) 


where M , is the molecular weight of species “z”, nr is the total number of reactions, v jt is the 

stoichiometric coefficient, K f . and K b . are the forward and backward reaction rates of the j th 

reaction, and ns is the total number of species. r' ] and r" are the power dependence of the 

reactant and product, respectively. The reaction rate, in general, is formulated in the empirical 
Arrhenius form as: 


K fj= A , TJ * 


B i„~ E a,jlRJ 


(E.3) 


where Ay, B,. E a j are empirical constants for j-th reaction. The input format is shown in the 
following table, where E/R is E a jlR u , STCOEF is v jt , and STCOEG is r ji . 
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Format 
(2 or 3 
lines* NR) 

REACTION: Species names, i = 1 — > ns (only 1 line as a title) 
j, A, B, E/R, ITHIRD, IGLOB, 

(STOCEF (/,;), / = 1 -> ns), 

(STOCEG (/,_/), i = 1 — » /z.v) * needed if IGLOB = 2 

Definition 

j 

reaction step counter 


A 

reaction rate leading constant 


B 

reaction rate temperature exponent 


E/R 

reaction rate activation energy constant 


ITHIRD 

third-body reaction indicator 
= 0: deactivated 

= n: for using the n-th species as third-body 

= 999: for global (every species) third-body 


IGLOB 

global reaction model indicator 
= 0: elementary reactions with the rate of backward 

reaction is calculated from equilibrium constant and 
forward reaction rate 

= 1: one-way reaction (either forward or backward 
reaction controlled by the sign of STCOEF); need 
only one input line of STCOEF 
= 2: one-way reaction with power dependency; need input 
line for STCOEF and STCOEG 


STCOEF 

stoichiometric coefficients of elementary reactions (negative 
signs apply to reactants and positive signs are for the 
products) 


STCOEG 

power dependency coefficients (negative signs apply to 
reactants and positive signs are for the products) 


For a sample reaction model: 


13 


REACTION :H20, 

02, 

CO, 

C02 , C7H8, 

H2 , 

o, 

H, 

OH, 

N2 

1, 4.4963E9, 

1 — * 

o 

o 

ro 

679E4, 

o, 

2, 





o., 

-3.5, 

7., 

o., 

-1. , 

4., 

o., 

o., 

o., 

0 . 

o., 

-1.0, 

o., 

0., -0.5, 

o., 

o., 

o., 

o., 

0 . 


2, 1.7000E13, 0.00,24070., 0, 0, 
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0., -1., 0., 0., 

0. 

r 

i., 

0., 

0., 

2., 

0 

3, 

2.1900E13, 0.00, 2590., 

0, 

o, 







1., 0., 0., 0., 

0. 

t ~ 

i., 

0., 

1., 

-1. , 

0 

4, 

6 . 0230E12 , 0.00, 550., 

0, 

0, 







1., 0., 0., 0., 

0. 

r 

o., 

1. , 

0., 

-2., 

0 

5, 

1 . 8000E10 , 1.00, 4480., 

0, 

o, 







o 

o 

o 

o 

0. 

t ~ 

i., 

-1. , 

1., 

1. , 

0 

6, 

1.2200E17, -0.91, 8369., 

0, 

0, 







0., -1., 0., 0., 

0. 

r 

0., 

1. , 

-1., 

1. , 

0 

7, 

4 . 0000E12 , 0.00,10030., 

0, 

o, 







0., 0., -1., 1., 

0. 

r 

0., 

0., 

1., 

-1. , 

0 

8, 

3.0000E12, 0.00,2. 50E4 , 

0, 

o. 







0., -1., -1., 1., 

0. 

r 

0., 

1. , 

0., 

0., 

0 

9, 

1 . 0000E16, 0.00, 0., 

999, 

o. 







o 

o 

o 

o 

0. 

r 

0., 

-1. , 

-1., 

1. , 

0 

o 
\ — 1 

2 . 5500E18 , -1.00,59390., 

999, 

o. 







0., 1., 0., 0., 

0. 

r 

0., 

-2., 

0., 

0., 

0 

11, 

5.0000E15, 0.00, 0., 

999, 

o. 







o 

o 

o 

o 

0. 

r 

1., 

0., 

-2., 

0., 

0 

12, 

8.4000E21, -2.00, 0., 

999, 

o. 







O 

O 

O 

0. 

r 

0., 

0., 

-1., 

-1. , 

0 

13, 

6.0000E13, 0.00, 0., 

999, 

o. 







0., 0., -1., 1., 

0. 

r 

0., 

-1. , 

0., 

0., 

0 


Reaction Step 

A 

B 

E a /R u 

Power Dependencies 
(r,j) 

C 7 H 8 + 3.502 7CO + 4H 2 

4.4963E9 

1 

2.679E4 

[C 7 H 8 ] 0 ' 5 [0 2 ]' 

H 2 + 0 2 <-> OH + OH 

1.7E13 

0 

2.407E4 

— 

OH + H 2 H 2 0 + H 

2.19E13 

0 

2.59E3 

— 

OH + OH 0 + H 2 0 

6.023E12 

0 

5.5E2 

— 

O + H 2 <h> H + OH 

1.8E10 

1.0 

4.48E3 

— 

H + 0 2 O + OH 

1.22E17 

-0.91 

8.369E3 

— 

H + 0+ M ^ OH + M 

1.0E16 

0 

0 

— 

o + o + m<->o 2 + m 

2.55E18 

-1.0 

5.939E4 

— 

h + h + m^h 2 + m 

5.0E15 

0 

0 

— 
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H + OH + Mo H 2 0 + M 

8.4E21 

-2.0 

0 

— 

CO + OH H + C0 2 

4.0E12 

0 

4.03E3 

— 

CO + O 2 0 CO 2 + 0 

3.0E12 

0 

2.5E4 

— 

CO + O + M C0 2 + M 

6.0E13 

0 

0 

— 


Note: M represents the third body, which includes all species in the above reaction model. If the 
power dependency is not specified, then the stoichiometric coefficient will be used. 
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